03-varcomp/VarianceComponentModels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\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",
    "**Lange Symposium**\n",
    "\n",
    "**Juhyun Kim, juhkim111@ucla.edu**\n",
    "\n",
    "**Department of Biostatistics, UCLA**\n",
    "\n",
    "**Feb 22, 2020**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Machine information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VarianceComponentModels.jl\n",
    "\n",
    "[`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",
    "$$\\text{vec}(Y) \\sim \\text{Normal}(XB, \\Sigma_1 \\otimes V_1 + \\cdots + \\Sigma_m \\otimes V_m)$$\n",
    "\n",
    "where $\\otimes$ is the [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product). \n",
    "\n",
    "In this model, **data** is represented by \n",
    "\n",
    "* `Y`: `n x d` response matrix \n",
    "* `X`: `n x p` covariate matrix \n",
    "* `V=(V1,...,Vm)`: a tuple of `m` `n x n` covariance matrices\n",
    "\n",
    "and **parameters** are\n",
    "\n",
    "* `B`: `p x d` mean parameter matrix\n",
    "* `Σ=(Σ1,...,Σm)`: a tuple of `m` `d x d` variance components. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Package Features \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)\n",
    "* [Heritability analysis](https://openmendel.github.io/VarianceComponentModels.jl/latest/man/heritability/#Heritability-Analysis-1) in genetics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heritability Analysis\n",
    "\n",
    "Variance component estimation can be used to estimate heritability of a quantitative trait. "
   ]
  },
  {
   "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. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SnpArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = normpath(SnpArrays.datadir())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`EUR_subset.bed`, `EUR_subset.bim`, and `EUR_subset.fam` is a set of Plink files in binary format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "maf_EUR = maf(EUR_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of minor allele frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## generate a histogram of MAF\n",
    "# using Plots, PyPlot\n",
    "# gr(size=(600, 500), html_output_format=:png)\n",
    "# hist_maf = histogram(maf_EUR, xlab = \"Minor Allele Frequency (MAF)\", \n",
    "#                    ylab = \"Number of SNPs\", label=\"\")\n",
    "# png(hist_maf, \"hist_MAF.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image info](./hist_MAF.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that about 29% of SNPs are rare variants (MAF < 0.05). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## GRM using SNPs with maf > 0.01 (default) \n",
    "Φgrm = grm(EUR_subset; method = :GRM) # classical genetic relationship matrix\n",
    "# Φ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 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`. \n",
    "\n",
    "The data should be available in `pheno.txt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## simulate `pheno.txt` \n",
    "# using LinearAlgebra, DelimitedFiles\n",
    "# Random.seed!(1234)\n",
    "# Ω = 0.1 * Φgrm + 0.9 * Matrix(1.0*I, nobs, nobs)\n",
    "# Ωchol = cholesky(Symmetric(Ω))\n",
    "# y = ones(nobs) + Ωchol.L * randn(nobs)\n",
    "# writedlm(\"pheno.txt\", y)"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DelimitedFiles \n",
    "pheno = readdlm(\"pheno.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram of phenotype values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## generate histogram of phenotype values\n",
    "#hist_pheno = histogram(pheno, xlab=\"Phenotype\", ylab=\"Frequency\", label=\"\")\n",
    "#png(hist_pheno, \"hist_pheno.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image info](./hist_pheno.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-processing data for heritability analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prepare variance component model fitting, we form an instance of VarianceComponentVariate. The two covariance matrices are $(2\\Phi, I)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VarianceComponentModels, LinearAlgebra\n",
    "# no. of observations \n",
    "nobs = size(pheno, 1)\n",
    "\n",
    "# form data as VarianceComponentVariate\n",
    "X = ones(nobs)\n",
    "EURdata = VarianceComponentVariate(pheno, X, (2Φgrm, Matrix(1.0I, nobs, nobs)))\n",
    "fieldnames(typeof(EURdata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EURdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Heritability 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 variane component estimation procedures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use Fisher scoring algorithm to fit variance component model for our trait. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ2a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Environmental/non-genetic variance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ2e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# heritability and its standard error from single trait analysis\n",
    "h, hse = heritability(trait_model.Σ, Σcov)\n",
    "\n",
    "\n",
    "[h[1], hse[1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also run MM algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": [
    "Heritability and its standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h, hse = heritability(trait_model.Σ, Σcov)\n",
    "[h[1], hse[1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate trait analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the joint analysis of multiple traits, go to [`VarianceComponentModels` documentation](https://openmendel.github.io/VarianceComponentModels.jl/latest/man/heritability/). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Repeat the above analysis computing the empirical kinship matrix using the method of moment method or the robust method. See [SnpArrays.jl documentation](https://openmendel.github.io/SnpArrays.jl/latest/#Genetic-relationship-matrix-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Φgrm = grm(EUR_subset; method = )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing SNP association using maximum likelihoods of variance component models\n",
    "credit: [Heritability 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 particular SNP has an effect on a given phenotype after accounting for relatedness among individuals. Here we fit variance component model with a single SNP *s* 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",
    "where \n",
    "\n",
    "* $\\mathbf{y}$: phenotype \n",
    "\n",
    "\n",
    "\n",
    "and \n",
    "\n",
    "\n",
    "* Fixed effects:\n",
    "    * $\\mathbf{X}$: matrix of covariates including intercept\n",
    "    * $\\beta$: vector of covariate effects, including intercept\n",
    "    * $\\mathbf{G}_s$: genotype of SNP *s*\n",
    "    * $\\gamma$: (scalar) 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 SNP *s* is associated with phenotype, we fit two models. First consider the model without SNP *s* 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 SNP *s* as fixed effects (1). Then we can compare the log likelihood to see if there is improvement in the model fit with inclusion of the SNP of interest. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the null model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VarianceComponentModels\n",
    "# null data model has two variance components but no SNP fixed effects\n",
    "# form data as VarianceComponentVariate matrix \n",
    "X = ones(nobs)\n",
    "nulldata = VarianceComponentVariate(pheno, X, (2Φgrm, Matrix(1.0I, nobs, nobs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nullmodel = VarianceComponentModel(nulldata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nulllogl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model mean effects (a grand mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nullmodel.B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model additive genetic variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nullmodel.Σ[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The null model environmental variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nullmodel.Σ[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit the alternative model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snp_vec = convert(Vector{Float64}, EUR_subset[:, 10]) \n",
    "Xalt = [ones(nobs) snp_vec]\n",
    "altdata = VarianceComponentVariate(pheno, Xalt, (2Φgrm, Matrix(1.0I, nobs, nobs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "altmodel = VarianceComponentModel(altdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@time altlogl, altmodel, = fit_mle!(altmodel, altdata; algo=:FS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model log-likelihood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "altlogl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model mean effects (a grand mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "altmodel.B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model additive genetic variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "altmodel.Σ[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The alternative model environmental variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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 1.786 (distributed chi-squared), with one degree of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Distributions\n",
    "LRT = 2(altlogl - nulllogl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The associated p-value: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pval = ccdf(Chisq(1), LRT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We see that adding this SNP as a covariate to the model does not fit significantly better than the null model. In other words, the SNP does 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.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.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
}