python_stats_intro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "plt.style.use('seaborn-talk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to pandas DataFrames: manipulating and exploring data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas is a Python package for easy to use data structures and analysis tools. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IQ and brain size data\n",
    "data = pd.read_csv('brain_size.csv', sep=';', na_values='.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DataFrames make it easy to subselect some of the data based on variable names. \n",
    "females = data[data['Gender']=='Female']\n",
    "females.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And to get some descriptive statistics out\n",
    "print('Mean female IQ: ')\n",
    "print( females['VIQ'].mean() )\n",
    "\n",
    "print('Standard devation: ')\n",
    "print(females['VIQ'].std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.groupby('Gender').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can easily create DataFrames from Numpy arrays:\n",
    "random_data = np.random.rand(20,2)\n",
    "data2 = pd.DataFrame(random_data, columns=['Height', 'Weight'])\n",
    "data2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And add other columns to your DataFrame, with a different datatype\n",
    "grades = ['A', 'B'] * 10\n",
    "data2['Grade'] = grades"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2.tail(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can change the index to your liking, for example here I'll reverse the index order\n",
    "data2.index = np.arange(20)[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# iloc gives the location (like a numpy array)\n",
    "data2.iloc[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2.loc[18]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting using pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pandas.tools import plotting\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The scatter matrix gives you a nice way to explore relations in your data (diagonals show histograms)\n",
    "plotting.scatter_matrix(data[['Weight', 'Height', 'MRI_Count', 'PIQ', 'FSIQ', 'VIQ']])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis testing using scipy\n",
    "\n",
    "scipy.stats is the go-to stats package for Python. It contains a large number of probability distributions, as well as a growing library of statistical functions and tests. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is the mean verbal IQ different from 100? the 1-sample t-test\n",
    "stats.ttest_1samp(data['VIQ'], 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Are the IQs of males and females different from each other? \n",
    "female_viq = data[data['Gender'] == 'Female']['VIQ']\n",
    "male_viq = data[data['Gender'] == 'Male']['VIQ']\n",
    "stats.ttest_ind(female_viq, male_viq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paired or dependent samples t-test for non-independent variables (here we compare different measures \n",
    "# of IQ within the same individuals)\n",
    "stats.ttest_rel(data['FSIQ'], data['PIQ'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The dependent samples t-test assumes normally distributed data. Scipy also has the standard non-parametric tests:\n",
    "stats.wilcoxon(data['FSIQ'], data['PIQ'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scipy.stats contains a wealth of probability distributions. \n",
    "\n",
    "from scipy.stats import invweibull\n",
    "fig, ax = plt.subplots(1,1)\n",
    "c = 10.6\n",
    "\n",
    "x = np.linspace(invweibull.ppf(0.01, c), invweibull.ppf(0.99, c), 100)\n",
    "ax.plot(x, invweibull.pdf(x, c),'k', lw=5, alpha=0.6, label='invweibull pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Wilcoxon sign test signed rank test is a close sibling of the dependent samples t-test.  Because the dependent samples t-tests analyzes if the average difference of two repeated measures is zero; it requires metric (interval or ratio) and normally distributed data; the Wilcoxon sign test uses ranked or ordinal data.  Thus it is a common alternative to the dependent samples t-test when its assumptions are not met."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statsmodels: \n",
    "## Linear models, ANOVA etc "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statsmodels is a module for the estimation of many different statistical models, conducting statistical tests and statistical data exploration. It is designed to work easily with Pandas dataframes and has an R-like syntax for defining the models. It has scipy.stats as a dependency and is complementatry to it rather than a substitute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's go back to our IQ data set:\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.formula.api import ols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ordinary least squares: \n",
    "# In statsmodels, we define the model using a formula like in R. Here, we try to predict IQ from Height and log Weight\n",
    "model = ols(\"FSIQ ~ MRI_Count\", data).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Similarly you can use the log of the MRI count like this:\n",
    "model = ols(\"FSIQ ~ Weight + np.log(MRI_Count)\", data).fit()\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical data: comparing groups or multiple categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Statsmodels can automatically infer a categorical variable. If you want to interpret a column of integers\n",
    "# you can force it to be categorical using C(). \n",
    "model = ols(\"VIQ ~ C(Gender)\", data).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in the iris data set (petal and sepal sizes of 3 different types of Iris flower)\n",
    "iris_data = pd.read_csv('iris.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "categories = pd.Categorical(iris_data['name'])\n",
    "categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.scatter_matrix(iris_data[['sepal_length', 'sepal_width', 'petal_length', 'petal_width']], c=categories.labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ols('sepal_width ~ name + petal_length + sepal_length', iris_data).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can do post-hoc testing for specific differences between our multiple predictors. \n",
    "# You can formulate contrasts with a list. Here we check the difference between versicolor and virginica:\n",
    "model.f_test([0, 1, -1, 0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing for interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ols('sepal_width ~ name + petal_length * petal_width', iris_data).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of variance (ANOVA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "table = sm.stats.anova_lm(model, typ=2)\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (By the way) Seaborn: easy visualisation and simple statistical fitting from pandas dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(iris_data, vars=['sepal_length','sepal_width','petal_length','petal_width'], kind='reg', diag_kind='kde')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(iris_data, vars=['sepal_length','sepal_width','petal_length','petal_width'], kind='reg', hue='name')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# rpy2 : Using R in python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation\n",
    "Installing on Ubuntu is super easy: install R with apt-get and pip install rpy2 \n",
    "\n",
    "For Mac I needed to install R through homebrew instead of the normal way (downloading from CRAN) in order for rpy2 to be able to find R. However, the homebrew version of R doesn't support X11, which is used for graphics. \n",
    "\n",
    "A workaround: \n",
    "- install X11 through homebrew:\n",
    "\n",
    "```\n",
    "brew cask install xquartz\n",
    "```\n",
    "- Install R with X11 support through this user's repo:\n",
    "\n",
    "```\n",
    "brew tap randy3k/r\n",
    "brew install r-x11\n",
    "```\n",
    "- Install rpy2 through pip:\n",
    "\n",
    "```\n",
    "pip install rpy2 \n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating R code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rpy2.robjects as robjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# robjects is your communication channel between Python and R. Basically, all rpy2 does is interface with an R\n",
    "# workspace and your python workspace. \n",
    "#help(robjects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can use the square brackets [] to get items from your R object, like the method __getitem__()\n",
    "pi = robjects.r['pi']\n",
    "pi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The result is an R vector, not a scalar. You can index it using the normal python way.\n",
    "pi[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# However the object is also callable. Any command that you pass through robjects.r('') will be interpreted as R code\n",
    "# For example:\n",
    "pi = robjects.r('pi')\n",
    "pi[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we first define a function in r, and then call it with the argument 3, which just gives us a float vector.\n",
    "robjects.r('''\n",
    "        f <- function(r, verbose=FALSE) {\n",
    "            if (verbose) {\n",
    "                cat(\"I am calling f().\\n\")\n",
    "            }\n",
    "            2 * pi * r\n",
    "        }\n",
    "        f(3)\n",
    "        ''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# however, in the r workspace, that function still exists, and we can make it accessible using our [] getitem method:\n",
    "r_f = robjects.r['f']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now use this as a regular python function:\n",
    "r_f(4)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "r_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# And we can see its R representation:\n",
    "print(r_f.r_repr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# watch out: because R represents everything as vectors, rather than scalars, we need to be careful.\n",
    "# The following code doesn't do quite what you'd expect:\n",
    "robjects.r['pi'] + 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# R String vectors: \n",
    "res = robjects.StrVector(['abc', 'def'])\n",
    "print(res.r_repr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Int vectors:\n",
    "res = robjects.IntVector([1, 2, 3])\n",
    "print(res.r_repr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Float vectors\n",
    "res = robjects.FloatVector([1.1, 2.2, 3.3])\n",
    "print(res.r_repr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# R matrixes and arrays are just vectors with a dim attribute.\n",
    "# The easiest way to create such objects is to do it through R functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = robjects.FloatVector([1.1, 2.2, 3.3, 4.4, 5.5, 6.6])\n",
    "m = robjects.r['matrix'](v, nrow = 2)\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calling R functions: \n",
    "rsum = robjects.r['sum']\n",
    "rsum(robjects.IntVector([1,2,3]))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting a function from R also allows you to use all of its keywords\n",
    "# Note: Arguments are now given in Python style, not R style. For example, a boolean will be 'True', not 'TRUE'\n",
    "rsort = robjects.r['sort']\n",
    "res = rsort(robjects.IntVector([1,2,3]), decreasing=True)\n",
    "print(res.r_repr())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphics and plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rpy2.robjects as robjects\n",
    "\n",
    "r = robjects.r\n",
    "\n",
    "x = robjects.IntVector(range(10))\n",
    "y = r.rnorm(10)\n",
    "r.X11()\n",
    "\n",
    "r.layout(r.matrix(robjects.IntVector([1,2,3,2]), nrow=2, ncol=2))\n",
    "r.plot(r.runif(10), y, xlab=\"runif\", ylab=\"foo/bar\", col=\"red\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2 import robjects\n",
    "from rpy2.robjects import Formula, Environment\n",
    "from rpy2.robjects.vectors import IntVector, FloatVector\n",
    "from rpy2.robjects.lib import grid\n",
    "from rpy2.robjects.packages import importr, data\n",
    "from rpy2.rinterface import RRuntimeError\n",
    "import warnings\n",
    "\n",
    "# The R 'print' function\n",
    "rprint = robjects.globalenv.get(\"print\")\n",
    "stats = importr('stats')\n",
    "grdevices = importr('grDevices')\n",
    "base = importr('base')\n",
    "datasets = importr('datasets')\n",
    "\n",
    "grid.activate()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can access ggplot like this:\n",
    "import math, datetime\n",
    "import rpy2.robjects.lib.ggplot2 as ggplot2\n",
    "import rpy2.robjects as ro\n",
    "from rpy2.robjects.packages import importr\n",
    "base = importr('base')\n",
    "\n",
    "mtcars = data(datasets).fetch('mtcars')['mtcars']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard R dataset containing fuel consumption and 10 aspects of automobile design and performance for 32 cars\n",
    "mtcars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can use ggplot to plot the miles per gallon (mpg) as a function of the weight of a car, separating out categories\n",
    "# of cars by their number of cylinders. \n",
    "pp = ggplot2.ggplot(mtcars) + \\\n",
    "     ggplot2.aes_string(x='wt', y='mpg', col='factor(cyl)') + \\\n",
    "     ggplot2.geom_point() + \\\n",
    "     ggplot2.geom_smooth(ggplot2.aes_string(group = 'cyl'),\n",
    "                         method = 'lm')\n",
    "pp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear models\n",
    "\n",
    "In R, specifying a linear model is very straightforward. Here's an example \n",
    "R code, where we want to predict the weight of a treatment and control group:\n",
    "\n",
    "```\n",
    "ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)\n",
    "trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)\n",
    "group <- gl(2, 10, 20, labels = c(\"Ctl\",\"Trt\"))\n",
    "weight <- c(ctl, trt)\n",
    "\n",
    "anova(lm.D9 <- lm(weight ~ group))\n",
    "\n",
    "summary(lm.D90 <- lm(weight ~ group - 1))# omitting intercept\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can do the same in rpy2 like this:\n",
    "\n",
    "from rpy2.robjects import FloatVector\n",
    "from rpy2.robjects.packages import importr\n",
    "\n",
    "stats = importr('stats')\n",
    "base = importr('base')\n",
    "\n",
    "ctl = FloatVector([4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14])\n",
    "trt = FloatVector([4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69])\n",
    "group = base.gl(2, 10, 20, labels = [\"Ctl\",\"Trt\"])\n",
    "weight = ctl + trt\n",
    "\n",
    "robjects.globalenv[\"weight\"] = weight\n",
    "robjects.globalenv[\"group\"] = group\n",
    "lm_D9 = stats.lm(\"weight ~ group\")\n",
    "print(stats.anova(lm_D9))\n",
    "\n",
    "# omitting the intercept\n",
    "lm_D90 = stats.lm(\"weight ~ group - 1\")\n",
    "print(base.summary(lm_D90))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How do we inspect the results? \n",
    "print(lm_D9.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can extract the R way using rx: \n",
    "lm_D9.rx2('coefficients')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or the Python way\n",
    "lm_D9[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal component analysis in rpy2 \n",
    "\n",
    "The R code for PCA on some random data is: \n",
    "```\n",
    "m <- matrix(rnorm(100), ncol=5)\n",
    "pca <- princomp(m)\n",
    "plot(pca, main=\"Eigen values\")\n",
    "biplot(pca, main=\"biplot\")\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In rpy2 this looks pretty similar\n",
    "import rpy2.robjects as robjects\n",
    "\n",
    "\n",
    "r = robjects.r\n",
    "r.x11()\n",
    "m = r.matrix(r.rnorm(100), ncol=5)\n",
    "pca = r.princomp(m)\n",
    "r.plot(pca, main=\"Eigen values\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r.x11()\n",
    "r.biplot(pca, main=\"biplot\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating multiple subplots\n",
    "\n",
    "from rpy2.robjects.packages import importr\n",
    "graphics = importr('graphics')\n",
    "grdevices = importr('grDevices')\n",
    "base = importr('base')\n",
    "stats = importr('stats')\n",
    "\n",
    "import array\n",
    "\n",
    "x = array.array('i', range(10))\n",
    "y = stats.rnorm(10)\n",
    "\n",
    "grdevices.X11()\n",
    "\n",
    "graphics.par(mfrow = array.array('i', [2,2]))\n",
    "graphics.plot(x, y, ylab = \"foo/bar\", col = \"red\")\n",
    "\n",
    "kwargs = {'ylab':\"foo/bar\", 'type':\"b\", 'col':\"blue\", 'log':\"x\"}\n",
    "graphics.plot(x, y, **kwargs)\n",
    "\n",
    "\n",
    "m = base.matrix(stats.rnorm(100), ncol=5)\n",
    "pca = stats.princomp(m)\n",
    "graphics.plot(pca, main=\"Eigen values\")\n",
    "stats.biplot(pca, main=\"biplot\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare rpy2 and scipy.stats: the wilcoxon test with very few data points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poplar_data = pd.read_csv('poplar.csv')\n",
    "poplar_data = poplar_data.drop('Unnamed: 0', axis=1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poplar_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only take first 7 observations: \n",
    "poplar_data = poplar_data.loc[:7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The scipy wilcoxon test uses a normal approximation to compute the p value from the test statistic. \n",
    "# This works for sufficient sample sizes, but with small sample sizes it gives problems and has no alternative. \n",
    "import scipy.stats as stats\n",
    "stats.wilcoxon(np.array(poplar_data['August']),np.array(poplar_data['November']), correction=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The equivalent function in R has the option to compute the exact p value instead of using the normal approximation\n",
    "# It will default to this option if the sample size is smaller than 50. \n",
    "# Using rpy2: \n",
    "from rpy2 import robjects\n",
    "from rpy2.robjects import FloatVector\n",
    "\n",
    "wilcox_r = robjects.r['wilcox.test']\n",
    "\n",
    "august = FloatVector(poplar_data['August'])\n",
    "november = FloatVector(poplar_data['November'])\n",
    "wilcox_r = robjects.r['wilcox.test']\n",
    "\n",
    "result = wilcox_r(august, november, paired=True, exact=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The interface with rpy2 can feel a bit clunky. If you work on a project with rpy2 you might want to \n",
    "# wrap your rpy2 commands in functions:\n",
    "\n",
    "def wilcoxon_r(x, y, **kwargs):\n",
    "    rx = FloatVector(x)\n",
    "    ry = FloatVector(y)\n",
    "    wilcox_r = robjects.r['wilcox.test']\n",
    "    result = wilcox_r(rx, ry, **kwargs)\n",
    "    statistic = result[0][0] # TODO: index by name instead\n",
    "    p_value = result[2][0]\n",
    "    return statistic, p_value\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This way you can use R functionality without having to think about these interfaces\n",
    "wilcoxon_r(poplar_data['August'], poplar_data['November'], paired=True, exact=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}