notebooks/VCM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VCM\n",
    "---\n",
    "\n",
    "Here, we have provided a sample example of using the variance component model.\n",
    "\n",
    "In this example we will use the [iris dataset](https://en.wikipedia.org/wiki/Iris_flower_data_set) as an example. We aim to investigate the proportion of variance in the petal length that can be explained by the observed inverted euclidean distance in petal width after regressing out sepal length and sepal width as nuisance variables.\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial imports\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from sklearn.metrics import pairwise_distances\n",
    "from statsmodels.stats.correlation_tools import cov_nearest\n",
    "\n",
    "# import the vcm code\n",
    "codes_directory = os.path.abspath('../codes')\n",
    "if codes_directory not in sys.path:\n",
    "    sys.path.append(codes_directory)\n",
    "import vcm\n",
    "\n",
    "# Note: you can import vcm directly if you copy the vcm.py code to your directory, this import method is only because we are accessing vcm.py from a different directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Variance ratio in petal_length explained by petal width: 0.91\n"
     ]
    }
   ],
   "source": [
    "# First, load the dataset:\n",
    "data = sns.load_dataset('iris')\n",
    "\n",
    "# variable of interest (trait values): petal_length\n",
    "y = np.array(data['petal_length'])\n",
    "\n",
    "# coveriates: sepal_width and sepal_length\n",
    "X = np.array(data[['sepal_length','sepal_width']])\n",
    "\n",
    "# Anatomical similarity: computed as inverted euclidean distance in petal_width\n",
    "euclidean_dist = pairwise_distances(np.array(data['petal_width']).reshape(-1, 1), metric = 'euclidean')\n",
    "K = cov_nearest(1 - (euclidean_dist/np.max(euclidean_dist)))\n",
    "\n",
    "# fit variance component model and report the ratio of overall variance explained by the similarity metric\n",
    "result = vcm.variance_component_model(y, X, K)\n",
    "\n",
    "print('--- Variance ratio in petal_length explained by petal width: {:.2f}'.format(result['m2']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}