{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# import some basic libraries\n",
    "import os\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# 3. ANOVA tables and post-hoc comparisons\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "<div class=\"alert alert-info\"><h4>Note</h4><p>ANOVAs and post-hoc tests are only available for `Lmer` models estimated using the `factors` argument of `model.fit()` and rely on implementations in R</p></div>\n",
    "\n",
    "In the previous tutorial where we looked at categorical predictors, behind the scenes `pymer4` was using the `factor` functionality in R. This means the output of `model.fit()` looks a lot like `summary()` in R applied to a model with categorical predictors. But what if we want to compute an F-test across *all levels* of our categorical predictor? \n",
    "\n",
    "`pymer4` makes this easy to do, and makes it easy to ensure Type III sums of squares infereces are valid. It also makes it easy to follow up omnibus tests with post-hoc pairwise comparisons. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## ANOVA tables and orthogonal contrasts\n",
    "\n",
    "Because ANOVA is just regression, `pymer4` can estimate ANOVA tables with F-results using the `.anova()` method on a fitted model. This will compute a Type-III SS table given the coding scheme provided when the model was initially fit. Based on the distribution of data across factor levels and the specific coding-scheme used, this may produce invalid Type-III SS computations. For this reason the `.anova()` method has a `force-orthogonal=True` argument that will reparameterize and refit the model using orthogonal polynomial contrasts prior to computing an ANOVA table.\n",
    "\n",
    "Here we first estimate a mode with dummy-coded categories and suppress the summary output of `.fit()`. Then we use `.anova()` to examine the F-test results. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SS Type III Analysis of Variance Table with Satterthwaite approximated degrees of freedom:\n",
      "(NOTE: Using original model contrasts, orthogonality not guaranteed)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SS</th>\n",
       "      <th>MS</th>\n",
       "      <th>NumDF</th>\n",
       "      <th>DenomDF</th>\n",
       "      <th>F-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>IV3</th>\n",
       "      <td>2359.778135</td>\n",
       "      <td>1179.889067</td>\n",
       "      <td>2</td>\n",
       "      <td>515.0</td>\n",
       "      <td>5.296284</td>\n",
       "      <td>0.005287</td>\n",
       "      <td>**</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              SS           MS  NumDF  DenomDF    F-stat     P-val Sig\n",
       "IV3  2359.778135  1179.889067      2    515.0  5.296284  0.005287  **"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# import model and sample data\n",
    "from pymer4.utils import get_resource_path\n",
    "from pymer4.models import Lmer\n",
    "\n",
    "# IV3 is a categorical predictors with 3 levels in the sample data\n",
    "df = pd.read_csv(os.path.join(get_resource_path(), 'sample_data.csv'))\n",
    "\n",
    "# # We're going to fit a multi-level regression using the \n",
    "# categorical predictor (IV3) which has 3 levels\n",
    "model = Lmer('DV ~ IV3 + (1|Group)', data=df)\n",
    "\n",
    "# Using dummy-coding; suppress summary output\n",
    "model.fit(factors={\n",
    "    'IV3': ['1.0', '0.5', '1.5']\n",
    "}, summarize=False)\n",
    "\n",
    "# Get ANOVA table\n",
    "model.anova()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Type III SS inferences will only be valid if data are fully balanced across levels or if contrasts between levels are orthogonally coded and sum to 0. Below we tell `pymer4` to respecify our contrasts to ensure this before estimating the ANOVA. `pymer4` also saves the last set of contrasts used priory to forcing orthogonality. \n",
    "\n",
    "Because the sample data is balanced across factor levels and there are not interaction terms, in this case orthogonal contrast coding doesn't change the results.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SS Type III Analysis of Variance Table with Satterthwaite approximated degrees of freedom:\n",
      "(NOTE: Model refit with orthogonal polynomial contrasts)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SS</th>\n",
       "      <th>MS</th>\n",
       "      <th>NumDF</th>\n",
       "      <th>DenomDF</th>\n",
       "      <th>F-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>IV3</th>\n",
       "      <td>2359.778135</td>\n",
       "      <td>1179.889067</td>\n",
       "      <td>2</td>\n",
       "      <td>515.000001</td>\n",
       "      <td>5.296284</td>\n",
       "      <td>0.005287</td>\n",
       "      <td>**</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              SS           MS  NumDF     DenomDF    F-stat     P-val Sig\n",
       "IV3  2359.778135  1179.889067      2  515.000001  5.296284  0.005287  **"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get ANOVA table, but this time force orthogonality \n",
    "# for valid SS III inferences\n",
    "# In this case the data are balanced so nothing changes\n",
    "model.anova(force_orthogonal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'IV3': ['0.5', '1.0', '1.5']}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Checkout current contrast scheme (for first contrast)\n",
    "# Notice how it's simply a linear contrast across levels\n",
    "model.factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'IV3': ['1.0', '0.5', '1.5']}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Checkout previous contrast scheme \n",
    "# which was a treatment contrast with 1.0\n",
    "# as the reference level\n",
    "model.factors_prev_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Marginal estimates and post-hoc comparisons\n",
    "-------------------------------------------\n",
    "`pymer4` leverages the `emmeans` package in order to compute marginal estimates (\"cell means\" in ANOVA lingo) and pair-wise comparisons of models that contain categorical terms and/or interactions. This can be performed by using the `.post_hoc()` method on fitted models. Let's see an example: \n",
    "\n",
    "First we'll quickly create a second categorical IV to demo with and estimate a 3x3 ANOVA to get main effects and the interaction.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SS Type III Analysis of Variance Table with Satterthwaite approximated degrees of freedom:\n",
      "(NOTE: Using original model contrasts, orthogonality not guaranteed)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SS</th>\n",
       "      <th>MS</th>\n",
       "      <th>NumDF</th>\n",
       "      <th>DenomDF</th>\n",
       "      <th>F-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>IV4</th>\n",
       "      <td>449.771051</td>\n",
       "      <td>224.885525</td>\n",
       "      <td>2</td>\n",
       "      <td>510.897775</td>\n",
       "      <td>1.006943</td>\n",
       "      <td>0.366058</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV3</th>\n",
       "      <td>2486.124318</td>\n",
       "      <td>1243.062159</td>\n",
       "      <td>2</td>\n",
       "      <td>508.993080</td>\n",
       "      <td>5.565910</td>\n",
       "      <td>0.004063</td>\n",
       "      <td>**</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV4:IV3</th>\n",
       "      <td>553.852530</td>\n",
       "      <td>138.463132</td>\n",
       "      <td>4</td>\n",
       "      <td>511.073624</td>\n",
       "      <td>0.619980</td>\n",
       "      <td>0.648444</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  SS           MS  NumDF     DenomDF    F-stat     P-val Sig\n",
       "IV4       449.771051   224.885525      2  510.897775  1.006943  0.366058    \n",
       "IV3      2486.124318  1243.062159      2  508.993080  5.565910  0.004063  **\n",
       "IV4:IV3   553.852530   138.463132      4  511.073624  0.619980  0.648444    "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fix the random number generator \n",
    "# for reproducibility\n",
    "import numpy as np\n",
    "np.random.seed(10)\n",
    "\n",
    "# Create a new categorical variable with 3 levels\n",
    "df = df.assign(IV4=np.random.choice(['1', '2', '3'], size=df.shape[0]))\n",
    "\n",
    "# Estimate model with orthogonal polynomial contrasts\n",
    "model = Lmer('DV ~ IV4*IV3 + (1|Group)', data=df)\n",
    "model.fit(factors={\n",
    "    'IV4': ['1', '2', '3'],\n",
    "    'IV3': ['1.0', '0.5', '1.5']},\n",
    "    ordered=True,\n",
    "    summarize=False\n",
    ")\n",
    "# Get ANOVA table\n",
    "# We can ignore the note in the output because\n",
    "# we manually specified polynomial contrasts\n",
    "model.anova()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Example 1\n",
    "\n",
    "Compare each level of IV3 to each other level of IV3, *within* each level of IV4. Use default Tukey HSD p-values.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-values adjusted by tukey method for family of 3 estimates\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>IV3</th>\n",
       "      <th>IV4</th>\n",
       "      <th>Estimate</th>\n",
       "      <th>2.5_ci</th>\n",
       "      <th>97.5_ci</th>\n",
       "      <th>SE</th>\n",
       "      <th>DF</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>42.554</td>\n",
       "      <td>33.778</td>\n",
       "      <td>51.330</td>\n",
       "      <td>4.398</td>\n",
       "      <td>68.140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.5</td>\n",
       "      <td>1</td>\n",
       "      <td>45.455</td>\n",
       "      <td>36.644</td>\n",
       "      <td>54.266</td>\n",
       "      <td>4.417</td>\n",
       "      <td>69.299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>40.904</td>\n",
       "      <td>32.196</td>\n",
       "      <td>49.612</td>\n",
       "      <td>4.361</td>\n",
       "      <td>65.943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2</td>\n",
       "      <td>42.092</td>\n",
       "      <td>33.301</td>\n",
       "      <td>50.882</td>\n",
       "      <td>4.406</td>\n",
       "      <td>68.609</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.5</td>\n",
       "      <td>2</td>\n",
       "      <td>41.495</td>\n",
       "      <td>32.829</td>\n",
       "      <td>50.161</td>\n",
       "      <td>4.339</td>\n",
       "      <td>64.626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>38.786</td>\n",
       "      <td>29.961</td>\n",
       "      <td>47.612</td>\n",
       "      <td>4.425</td>\n",
       "      <td>69.746</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>43.424</td>\n",
       "      <td>34.741</td>\n",
       "      <td>52.107</td>\n",
       "      <td>4.348</td>\n",
       "      <td>65.149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.5</td>\n",
       "      <td>3</td>\n",
       "      <td>46.008</td>\n",
       "      <td>37.261</td>\n",
       "      <td>54.755</td>\n",
       "      <td>4.383</td>\n",
       "      <td>67.208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>38.119</td>\n",
       "      <td>29.384</td>\n",
       "      <td>46.854</td>\n",
       "      <td>4.376</td>\n",
       "      <td>66.801</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   IV3 IV4  Estimate  2.5_ci  97.5_ci     SE      DF\n",
       "1  1.0   1    42.554  33.778   51.330  4.398  68.140\n",
       "2  0.5   1    45.455  36.644   54.266  4.417  69.299\n",
       "3  1.5   1    40.904  32.196   49.612  4.361  65.943\n",
       "4  1.0   2    42.092  33.301   50.882  4.406  68.609\n",
       "5  0.5   2    41.495  32.829   50.161  4.339  64.626\n",
       "6  1.5   2    38.786  29.961   47.612  4.425  69.746\n",
       "7  1.0   3    43.424  34.741   52.107  4.348  65.149\n",
       "8  0.5   3    46.008  37.261   54.755  4.383  67.208\n",
       "9  1.5   3    38.119  29.384   46.854  4.376  66.801"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute post-hoc tests\n",
    "marginal_estimates, comparisons = model.post_hoc(marginal_vars='IV3', grouping_vars='IV4')\n",
    "\n",
    "# \"Cell\" means of the ANOVA\n",
    "marginal_estimates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Contrast</th>\n",
       "      <th>IV4</th>\n",
       "      <th>Estimate</th>\n",
       "      <th>2.5_ci</th>\n",
       "      <th>97.5_ci</th>\n",
       "      <th>SE</th>\n",
       "      <th>DF</th>\n",
       "      <th>T-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>1</td>\n",
       "      <td>-2.901</td>\n",
       "      <td>-9.523</td>\n",
       "      <td>3.721</td>\n",
       "      <td>2.817</td>\n",
       "      <td>510.016</td>\n",
       "      <td>-1.030</td>\n",
       "      <td>0.558</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>1.650</td>\n",
       "      <td>-4.750</td>\n",
       "      <td>8.050</td>\n",
       "      <td>2.723</td>\n",
       "      <td>510.137</td>\n",
       "      <td>0.606</td>\n",
       "      <td>0.817</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>4.552</td>\n",
       "      <td>-1.951</td>\n",
       "      <td>11.054</td>\n",
       "      <td>2.766</td>\n",
       "      <td>510.267</td>\n",
       "      <td>1.645</td>\n",
       "      <td>0.228</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>2</td>\n",
       "      <td>0.596</td>\n",
       "      <td>-5.749</td>\n",
       "      <td>6.942</td>\n",
       "      <td>2.700</td>\n",
       "      <td>510.249</td>\n",
       "      <td>0.221</td>\n",
       "      <td>0.973</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>3.305</td>\n",
       "      <td>-3.387</td>\n",
       "      <td>9.998</td>\n",
       "      <td>2.847</td>\n",
       "      <td>510.883</td>\n",
       "      <td>1.161</td>\n",
       "      <td>0.477</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>2.709</td>\n",
       "      <td>-3.749</td>\n",
       "      <td>9.166</td>\n",
       "      <td>2.747</td>\n",
       "      <td>510.732</td>\n",
       "      <td>0.986</td>\n",
       "      <td>0.586</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>3</td>\n",
       "      <td>-2.584</td>\n",
       "      <td>-8.893</td>\n",
       "      <td>3.725</td>\n",
       "      <td>2.684</td>\n",
       "      <td>510.213</td>\n",
       "      <td>-0.963</td>\n",
       "      <td>0.601</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>5.305</td>\n",
       "      <td>-1.006</td>\n",
       "      <td>11.615</td>\n",
       "      <td>2.685</td>\n",
       "      <td>510.710</td>\n",
       "      <td>1.976</td>\n",
       "      <td>0.119</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>7.889</td>\n",
       "      <td>1.437</td>\n",
       "      <td>14.340</td>\n",
       "      <td>2.745</td>\n",
       "      <td>510.663</td>\n",
       "      <td>2.874</td>\n",
       "      <td>0.012</td>\n",
       "      <td>*</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Contrast IV4  Estimate  2.5_ci  97.5_ci     SE       DF  T-stat  P-val Sig\n",
       "1  1.0 - 0.5   1    -2.901  -9.523    3.721  2.817  510.016  -1.030  0.558    \n",
       "2  1.0 - 1.5   1     1.650  -4.750    8.050  2.723  510.137   0.606  0.817    \n",
       "3  0.5 - 1.5   1     4.552  -1.951   11.054  2.766  510.267   1.645  0.228    \n",
       "4  1.0 - 0.5   2     0.596  -5.749    6.942  2.700  510.249   0.221  0.973    \n",
       "5  1.0 - 1.5   2     3.305  -3.387    9.998  2.847  510.883   1.161  0.477    \n",
       "6  0.5 - 1.5   2     2.709  -3.749    9.166  2.747  510.732   0.986  0.586    \n",
       "7  1.0 - 0.5   3    -2.584  -8.893    3.725  2.684  510.213  -0.963  0.601    \n",
       "8  1.0 - 1.5   3     5.305  -1.006   11.615  2.685  510.710   1.976  0.119    \n",
       "9  0.5 - 1.5   3     7.889   1.437   14.340  2.745  510.663   2.874  0.012   *"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Pairwise comparisons\n",
    "comparisons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Example 2\n",
    "\n",
    "Compare each unique IV3,IV4 \"cell mean\" to every other IV3,IV4 \"cell mean\" and used FDR correction for multiple comparisons:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-values adjusted by fdr method for 36 comparisons\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Contrast</th>\n",
       "      <th>Estimate</th>\n",
       "      <th>2.5_ci</th>\n",
       "      <th>97.5_ci</th>\n",
       "      <th>SE</th>\n",
       "      <th>DF</th>\n",
       "      <th>T-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0,1 - 0.5,1</td>\n",
       "      <td>-2.901</td>\n",
       "      <td>-11.957</td>\n",
       "      <td>6.155</td>\n",
       "      <td>2.817</td>\n",
       "      <td>510.016</td>\n",
       "      <td>-1.030</td>\n",
       "      <td>0.535</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0,1 - 1.5,1</td>\n",
       "      <td>1.650</td>\n",
       "      <td>-7.102</td>\n",
       "      <td>10.403</td>\n",
       "      <td>2.723</td>\n",
       "      <td>510.137</td>\n",
       "      <td>0.606</td>\n",
       "      <td>0.726</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0,1 - 1.0,2</td>\n",
       "      <td>0.463</td>\n",
       "      <td>-8.657</td>\n",
       "      <td>9.582</td>\n",
       "      <td>2.837</td>\n",
       "      <td>511.103</td>\n",
       "      <td>0.163</td>\n",
       "      <td>0.871</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0,1 - 0.5,2</td>\n",
       "      <td>1.059</td>\n",
       "      <td>-7.649</td>\n",
       "      <td>9.766</td>\n",
       "      <td>2.709</td>\n",
       "      <td>510.435</td>\n",
       "      <td>0.391</td>\n",
       "      <td>0.835</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.0,1 - 1.5,2</td>\n",
       "      <td>3.768</td>\n",
       "      <td>-5.364</td>\n",
       "      <td>12.899</td>\n",
       "      <td>2.841</td>\n",
       "      <td>510.737</td>\n",
       "      <td>1.326</td>\n",
       "      <td>0.473</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Contrast  Estimate  2.5_ci  97.5_ci     SE       DF  T-stat  P-val Sig\n",
       "1  1.0,1 - 0.5,1    -2.901 -11.957    6.155  2.817  510.016  -1.030  0.535    \n",
       "2  1.0,1 - 1.5,1     1.650  -7.102   10.403  2.723  510.137   0.606  0.726    \n",
       "3  1.0,1 - 1.0,2     0.463  -8.657    9.582  2.837  511.103   0.163  0.871    \n",
       "4  1.0,1 - 0.5,2     1.059  -7.649    9.766  2.709  510.435   0.391  0.835    \n",
       "5  1.0,1 - 1.5,2     3.768  -5.364   12.899  2.841  510.737   1.326  0.473    "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute post-hoc tests\n",
    "marginal_estimates, comparisons = model.post_hoc(marginal_vars=['IV3', 'IV4'], p_adjust='fdr')\n",
    "\n",
    "# Pairwise comparisons, print the head because there are a lot\n",
    "# marginal_estimates are the same as the previous example\n",
    "comparisons.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Example 3\n",
    "\n",
    "For this example we'll estimate a more complicated ANOVA with 1 continuous IV and 2 categorical IVs with 3 levels each. This is the same model as before but with IV2 thrown into the mix. Now, pairwise comparisons reflect changes in the *slope* of the continuous IV (IV2) between levels of the categorical IVs (IV3 and IV4).\n",
    "\n",
    "First let's get the ANOVA table\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SS Type III Analysis of Variance Table with Satterthwaite approximated degrees of freedom:\n",
      "(NOTE: Using original model contrasts, orthogonality not guaranteed)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SS</th>\n",
       "      <th>MS</th>\n",
       "      <th>NumDF</th>\n",
       "      <th>DenomDF</th>\n",
       "      <th>F-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>IV2</th>\n",
       "      <td>46010.245471</td>\n",
       "      <td>46010.245471</td>\n",
       "      <td>1</td>\n",
       "      <td>535.763367</td>\n",
       "      <td>306.765451</td>\n",
       "      <td>1.220547e-54</td>\n",
       "      <td>***</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV3</th>\n",
       "      <td>726.318000</td>\n",
       "      <td>363.159000</td>\n",
       "      <td>2</td>\n",
       "      <td>500.573997</td>\n",
       "      <td>2.421301</td>\n",
       "      <td>8.984551e-02</td>\n",
       "      <td>.</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV4</th>\n",
       "      <td>143.379932</td>\n",
       "      <td>71.689966</td>\n",
       "      <td>2</td>\n",
       "      <td>502.297291</td>\n",
       "      <td>0.477981</td>\n",
       "      <td>6.203159e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV2:IV3</th>\n",
       "      <td>613.455876</td>\n",
       "      <td>306.727938</td>\n",
       "      <td>2</td>\n",
       "      <td>500.403443</td>\n",
       "      <td>2.045056</td>\n",
       "      <td>1.304528e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV2:IV4</th>\n",
       "      <td>4.914900</td>\n",
       "      <td>2.457450</td>\n",
       "      <td>2</td>\n",
       "      <td>502.300664</td>\n",
       "      <td>0.016385</td>\n",
       "      <td>9.837494e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV3:IV4</th>\n",
       "      <td>92.225327</td>\n",
       "      <td>23.056332</td>\n",
       "      <td>4</td>\n",
       "      <td>502.950771</td>\n",
       "      <td>0.153724</td>\n",
       "      <td>9.612985e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>IV2:IV3:IV4</th>\n",
       "      <td>368.085569</td>\n",
       "      <td>92.021392</td>\n",
       "      <td>4</td>\n",
       "      <td>503.354865</td>\n",
       "      <td>0.613537</td>\n",
       "      <td>6.530638e-01</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       SS            MS  NumDF     DenomDF      F-stat         P-val  Sig\n",
       "IV2          46010.245471  46010.245471      1  535.763367  306.765451  1.220547e-54  ***\n",
       "IV3            726.318000    363.159000      2  500.573997    2.421301  8.984551e-02    .\n",
       "IV4            143.379932     71.689966      2  502.297291    0.477981  6.203159e-01     \n",
       "IV2:IV3        613.455876    306.727938      2  500.403443    2.045056  1.304528e-01     \n",
       "IV2:IV4          4.914900      2.457450      2  502.300664    0.016385  9.837494e-01     \n",
       "IV3:IV4         92.225327     23.056332      4  502.950771    0.153724  9.612985e-01     \n",
       "IV2:IV3:IV4    368.085569     92.021392      4  503.354865    0.613537  6.530638e-01     "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = Lmer('DV ~ IV2*IV3*IV4 + (1|Group)', data=df)\n",
    "# Only need to polynomial contrasts for IV3 and IV4\n",
    "# because IV2 is continuous\n",
    "model.fit(factors={\n",
    "    'IV4': ['1', '2', '3'],\n",
    "    'IV3': ['1.0', '0.5', '1.5']},\n",
    "    ordered=True,\n",
    "    summarize=False\n",
    ")\n",
    "\n",
    "# Get ANOVA table\n",
    "model.anova()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Now we can compute the pairwise difference in slopes \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "Collapsed": "false",
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-values adjusted by bonf method for 3 comparisons\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Contrast</th>\n",
       "      <th>IV4</th>\n",
       "      <th>Estimate</th>\n",
       "      <th>2.5_ci</th>\n",
       "      <th>97.5_ci</th>\n",
       "      <th>SE</th>\n",
       "      <th>DF</th>\n",
       "      <th>T-stat</th>\n",
       "      <th>P-val</th>\n",
       "      <th>Sig</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.053</td>\n",
       "      <td>-0.254</td>\n",
       "      <td>0.147</td>\n",
       "      <td>0.084</td>\n",
       "      <td>502.345</td>\n",
       "      <td>-0.638</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.131</td>\n",
       "      <td>-0.313</td>\n",
       "      <td>0.050</td>\n",
       "      <td>0.076</td>\n",
       "      <td>502.494</td>\n",
       "      <td>-1.734</td>\n",
       "      <td>0.250</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.078</td>\n",
       "      <td>-0.278</td>\n",
       "      <td>0.122</td>\n",
       "      <td>0.083</td>\n",
       "      <td>502.821</td>\n",
       "      <td>-0.933</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>2</td>\n",
       "      <td>-0.038</td>\n",
       "      <td>-0.210</td>\n",
       "      <td>0.134</td>\n",
       "      <td>0.072</td>\n",
       "      <td>501.096</td>\n",
       "      <td>-0.526</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>0.002</td>\n",
       "      <td>-0.184</td>\n",
       "      <td>0.189</td>\n",
       "      <td>0.078</td>\n",
       "      <td>502.745</td>\n",
       "      <td>0.031</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>2</td>\n",
       "      <td>0.040</td>\n",
       "      <td>-0.142</td>\n",
       "      <td>0.222</td>\n",
       "      <td>0.076</td>\n",
       "      <td>502.836</td>\n",
       "      <td>0.530</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.0 - 0.5</td>\n",
       "      <td>3</td>\n",
       "      <td>-0.134</td>\n",
       "      <td>-0.329</td>\n",
       "      <td>0.061</td>\n",
       "      <td>0.081</td>\n",
       "      <td>502.956</td>\n",
       "      <td>-1.646</td>\n",
       "      <td>0.301</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.0 - 1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>-0.110</td>\n",
       "      <td>-0.302</td>\n",
       "      <td>0.083</td>\n",
       "      <td>0.080</td>\n",
       "      <td>502.109</td>\n",
       "      <td>-1.368</td>\n",
       "      <td>0.516</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.5 - 1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>0.024</td>\n",
       "      <td>-0.166</td>\n",
       "      <td>0.214</td>\n",
       "      <td>0.079</td>\n",
       "      <td>502.538</td>\n",
       "      <td>0.304</td>\n",
       "      <td>1.000</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Contrast IV4  Estimate  2.5_ci  97.5_ci     SE       DF  T-stat  P-val Sig\n",
       "1  1.0 - 0.5   1    -0.053  -0.254    0.147  0.084  502.345  -0.638  1.000    \n",
       "2  1.0 - 1.5   1    -0.131  -0.313    0.050  0.076  502.494  -1.734  0.250    \n",
       "3  0.5 - 1.5   1    -0.078  -0.278    0.122  0.083  502.821  -0.933  1.000    \n",
       "4  1.0 - 0.5   2    -0.038  -0.210    0.134  0.072  501.096  -0.526  1.000    \n",
       "5  1.0 - 1.5   2     0.002  -0.184    0.189  0.078  502.745   0.031  1.000    \n",
       "6  0.5 - 1.5   2     0.040  -0.142    0.222  0.076  502.836   0.530  1.000    \n",
       "7  1.0 - 0.5   3    -0.134  -0.329    0.061  0.081  502.956  -1.646  0.301    \n",
       "8  1.0 - 1.5   3    -0.110  -0.302    0.083  0.080  502.109  -1.368  0.516    \n",
       "9  0.5 - 1.5   3     0.024  -0.166    0.214  0.079  502.538   0.304  1.000    "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute post-hoc tests with bonferroni correction\n",
    "marginal_estimates, comparisons = model.post_hoc(marginal_vars='IV2',\n",
    "                                    grouping_vars=['IV3', 'IV4'],\n",
    "                                    p_adjust='bonf')\n",
    "\n",
    "# Pairwise comparisons\n",
    "comparisons"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}