Two-Way and Multi-way ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-way ANOVA\n",
    "For this first time in this course, we're going to talk about having multiple IVs that will predict your DV.  In the case of multi-way (or factorial) ANOVA all of the IVs should be categorical/groups and will predict our numerical DV.  We can also do a repeated measures ANOVA with additional IV(s) that are grouping variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "[Basic Two-way ANOVA](#basic)\n",
    "    <ul>\n",
    "    <li>[Check assumptions](#precheck) </li>\n",
    "    <li>[Main Effects](#mainanova)</li>\n",
    "    <li>[Interactions](#interact)</li>\n",
    "    <li>[Post-hoc Normality of Residuals](#residnorm)</li>\n",
    "    <li>[Post-hoc Pairwise t-tests](#pairwise)</li>\n",
    "    <li>[Effect Size](#rsquare)</li>\n",
    "    </ul>\n",
    "     <BR>\n",
    "[Repeated Measures](#repeat)\n",
    "    <ul>\n",
    "    <li>[Setup and Check Assumptions](#setup)</li>\n",
    "    <li>[The ANOVA](#repeatanova)</li>\n",
    "    <li>[Post-hoc Normality of Residuals](#residnormrep)</li>\n",
    "    <li>[Effect Size](#rsquarerep)</li>\n",
    "        </ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## load the libraries we'll need in the notebook\n",
    "\n",
    "library(tidyverse)\n",
    "library(magrittr)\n",
    "library(ggpubr) # containes line/dot plot for visualizing means\n",
    "library(DescTools) # contains levene's test function\n",
    "library(pwr) # for power analysis\n",
    "library(gridExtra) # for side-by-side plots\n",
    "library(rsq) # for r-squared functions\n",
    "library(afex) # for repeated measures ANOVA\n",
    "\n",
    "options(repr.plot.width=10, repr.plot.height=4) ## set options for plot size within the notebook -\n",
    "# this is only for jupyter notebooks, you can disregard this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"basic\"></a>\n",
    "## Basic Two-way ANOVA\n",
    "\n",
    "For the two-way ANOVA we are going to use the built-in dataset \"warpbreaks.\"  This dataset has three variables that describe warpbreaks on looms (it has something to do with yarn and weaving.  The unit of observation (row) are the looms.  The varibles are: <br>\n",
    "[,1]\tbreaks\t numeric\t The number of breaks <BR>\n",
    "[,2]\twool\t factor\t The type of wool (A or B)<BR>\n",
    "[,3]\ttension\t factor\t The level of tension (L, M, H)<BR>\n",
    "<BR>\n",
    "Let's start by looking at our variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data(warpbreaks)\n",
    "summary(warpbreaks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to look at grouped boxplots for both of our categorical variables.  We will use grid.arrange() from the package \"gridExtra\" so that we can place the plots side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#build first plot and save object\n",
    "wool_box <- warpbreaks %>% ggplot(aes(x = wool, y = breaks, fill = wool)) +\n",
    "        geom_boxplot()\n",
    "#build second plot and save object\n",
    "tens_box <- warpbreaks %>% ggplot(aes(x = tension, y = breaks, fill = tension)) +\n",
    "        geom_boxplot()\n",
    "\n",
    "#arrange saved plots side-by-side\n",
    "grid.arrange(wool_box, tens_box, ncol = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we'll also create grouped mean plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#first plot\n",
    "wool_plot <- warpbreaks %>% ggline(x = \"wool\", y = \"breaks\", \n",
    "       add = c(\"mean_se\", \"jitter\"),  add.params = list(color=\"wool\"),\n",
    "       ylab = \"Breaks\", xlab = \"Type of Wool\") \n",
    "\n",
    "#second plot\n",
    "tens_plot <- warpbreaks %>% ggline(x = \"tension\", y = \"breaks\", \n",
    "       add = c(\"mean_se\", \"jitter\"),  add.params = list(color=\"tension\"),\n",
    "       ylab = \"Breaks\", xlab = \"Tension Setting\") \n",
    "\n",
    "#arrange and print\n",
    "grid.arrange(wool_plot, tens_plot, ncol = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"precheck\"></a>\n",
    "### Check Assumptions\n",
    "- DV is numeric (interval or ratio) - yes\n",
    "- No extreme outliers - nothing looks too bad\n",
    "- Normality __*of residuals*__ - can't check this until after\n",
    "- Independence of Observations (random selection, different samples) - yes - they each come from different looms\n",
    "- Group sample sizes are approximately equal - exactly equal\n",
    "\n",
    "AND....\n",
    "- Homogeneity of Variance - Let's check this right now with Levene's test.\n",
    "\n",
    "Recall:\n",
    "\n",
    "$H_0:$ The variances in the groups are equal. <BR>\n",
    "$H_A:$ The variances in the groups are not equal.\n",
    "\n",
    "In this test, we sort of want to fail to reject null, because it's easier if our variances are equal and we don't need to make the adjustment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#LeveneTest(DV ~ IV, data = your data frame)\n",
    "\n",
    "LeveneTest(breaks ~ tension, data = warpbreaks)\n",
    "LeveneTest(breaks ~ wool, data = warpbreaks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"mainanova\"></a>\n",
    "### The Main Effects\n",
    "First we will run an ANOVA that looks at the main effects of each of the grouping variables.  This is similar to the one-way ANOVA, but instead of just separating the between variance of the one set of groups vs. the residuals, we have the between group variance for two sets of groups, and then the left over residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "summary(aov(breaks ~ wool + tension, data=warpbreaks))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean number of breaks differs by tension, but not by wool.  So we would say that tension, but not wool, influences breaks in warp on looms.\n",
    "\n",
    "Notice we get 2 F-values, this means that we have two sets of null/alternative hypotheses, one for each IV.\n",
    "\n",
    "For wool:<br>\n",
    "$H_0$: The means of breaks between all the wool groupings are equal.<br>\n",
    "$H_A$: The means of breaks between the wool groupings are not equal; there is at least one mean that is significantly different.<br>\n",
    "<br>\n",
    "For tension:<br>\n",
    "$H_0$: The means of breaks between the tension groups are equal.<br>\n",
    "$H_A$: The means of breaks between the tension groups are not equal; there is at least one mean that is significantly different.<br>\n",
    "\n",
    "Note the degrees of freedom.  There are 54 observations in the warpbreaks data.  The total degrees of freedom remains n-1.  wool gets 1 degree of freedom - $k_1 - 1$ where $k_1$ is the number of groups/levels of wool, tension gets 2 degrees of freedom ($k_2 - 1$ where $k_2 = 3$) and the remaining degrees of freedom are the residual degrees of freedom = $n - k_1 - k_2 - 1 = 50$\n",
    "\n",
    "<a id=\"interact\"></a>\n",
    "\n",
    "### Interaction Effect\n",
    "In some cases there is shared variance explained by the combination of both grouping variables.  This is called an interaction effect.  The effect of wool on breaks differs by the tension (or could be explained in the opposite way - the effect of tenion on breaks differs by the type of wool.  So there are the individual main effects - the variance that is explained exclusively by that variable, and the additional interaction effect that is the shared variance between the two grouping variables.\n",
    "\n",
    "Let's look at a graph in which we overlay both the means by wool and means by tension on the same graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## THIS OPTIONS IS ONLY FOR BINDER DISPLAY!!!\n",
    "options(repr.plot.width=5, repr.plot.height=4)\n",
    "\n",
    "#first plot\n",
    "warpbreaks %>% ggline(x = \"tension\", y = \"breaks\", color = \"wool\",\n",
    "       add = c(\"mean_se\", \"jitter\"),\n",
    "       palette = c(\"#00BF7D\", \"#FF61C9\"), size = 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The place in which the lines cross over each other is the place in which wool and tension interact.  When the tension is either low or high wool A has more breaks, however when there is medium tension, wool B has a higher mean number of breaks.  Adding in that interaction effect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "summary(aov(breaks ~ wool * tension, data=warpbreaks)) ## if you don't print summary you don't get full output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that tension remains significant, and in addition wool:tension interaction is significant.  Although wool is not a significant predictor of breaks on its own, the interplay of wool and tension is significant - wool acts differently depending on the tension used.\n",
    "\n",
    "Now that we have a third F-value, we have a third set of null/alternative hypotheses.  In addition to those we specified above for wool and tension separately, we now have a set for the interaction.\n",
    "\n",
    "For wool:tension (the interaction term) - <br>\n",
    "$H_0$: There is no interaction between wool and tension.<br>\n",
    "$H_A$: There is an interaction between wool and tension.<br>\n",
    "\n",
    "Here the degrees of freedom are different - the interaction gets some of the degrees of freedom.  There are 54 observations in the warpbreaks data.  The total degrees of freedom remains n-1 = 53.  \n",
    "\n",
    "wool gets 1 degree of freedom - $k_1 - 1$ where $k_1$ is the number of groups/levels of wool\n",
    "\n",
    "tension gets 2 degrees of freedom - $k_2 - 1$ where $k_2 = 3$ \n",
    "\n",
    "wool:tension gets 2 degrees of freedom $(k_1 - 1)(k_2 - 1) = (2-1)(3-1) = 2$\n",
    "\n",
    "and the remaining degrees of freedom are the residual degrees of freedom = $n - k_1 - k_2 - (k_1 - 1)(k_2 - 1) - 1 = 48$\n",
    "\n",
    "It's time for our post hoc analyses.\n",
    "\n",
    "<a id=\"residnorm\"></a>\n",
    "### Post-hoc check of the normality of the residuals\n",
    "As always, we need to check the normality of our residuals.  We will check this for our interaction ANOVA.  Remember that we inspect this graphically via a QQ plot, which is the second output when we plot a saved ANOVA object with the plot() function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## save interaction anova to object\n",
    "wb_aov <- aov(breaks ~ wool * tension, data=warpbreaks)\n",
    "\n",
    "## plot the post-hoc graphs\n",
    "plot(wb_aov)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residuals in the second plot follow closely along that reference diagonal, therefore the residuals are normally distributed.\n",
    "\n",
    "Now that we've dealt with our post-hoc assumption, we can move on to post-hoc analyses - pair-wise tests.\n",
    "\n",
    "<a id=\"pairwise\"></a>\n",
    "### Post-Hoc Pairwise Tests\n",
    "We can run our TukeyHSD() analysis on our saved aov object, the same way we did with our one-way ANOVA, however this time we will get 3 sets of pairwise tests - pairwise tests for just wool, pairwise tests for just tension, and pairwise tests for the combinations (interactions) of wool and tension.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "TukeyHSD(wb_aov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wool had only one pairwise comparison - A vs. B, since it only had two levels.  Tension has 3 pairwise comparisons for its three groups and only two are significant (M vs L and H vs L).  Then we have all the possible combinations of all the possible different levels of wool:tension.  The first 5 are significant.\n",
    "<a id=\"rsquare\"></a>\n",
    "### Effect Size - R-squared\n",
    "Finally, we want to know how substantively significant our result is.  In this situation we have 3 things that are explaining the variance in our DV(breaks) - 1) the main effect of wool, 2) the main effect of tension, 3) the interaction effect of wool:tension.  When we run these models the most typical output is the overall r-squared, which reflects the entire proportion of variance explained by all of the IVs (and interactions).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "wb_lm <- lm(breaks ~ wool * tension, data = warpbreaks)\n",
    "summary(wb_lm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our overall ANOVA has an r-squared of 0.3778, which means that 38% of the variance in breaks is explained by wool, tension, and the interaction of wool:tension.  If we want to know what the effect of each IV is individually, we can use eta-squared and partial eta-squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta2 <- EtaSq(wb_aov) ## give eta-squared the saved anova output\n",
    "eta2 # print the entire eta-squared output\n",
    "sum(eta2[1:3]) # sum the first column (eta.sq) to show that it equals the r-squared (unadjusted) above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see the component parts of the variance explained by each IV.  The first column, labeled eta.sq, portions out the portion of the IV over the total variance.  The partial eta-squared (eta.sq.part) is the SS of that IV over the SS of that IV plus just the SS of the residuals.  \n",
    "\n",
    "# $eta.sq_{wool} = \\frac{SS_{wool}}{SS_{total}}$   &nbsp;&nbsp;&nbsp;   but   &nbsp;&nbsp;  &nbsp; $eta.sq.part_{wool} = \\frac{SS_{wool}}{SS_{wool} + SS_{resid}}$\n",
    "\n",
    "This is the effect of the variable in isolation of the other predictors.  This value is not necessarily useful in a situation with a single analysis, but is useful in comparing the effect of a certain predictor across experiments where there are different IVs included in each.\n",
    "\n",
    "Note: Eta-squared is always a slightly biased (inflated) estimate of population variance explained, which is why you will see lm() output has both an r-squared and an \"adjusted r-squared\" reported."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"repeat\"></a>\n",
    "## Repeated Measures ANOVA, with additional IV\n",
    "Another form of two-way (or multi-way) ANOVA is to combine a repeated measure (time predictor) with a second, not time-related, categorical IV.  For this example we'll look at data from the built-in dataset \"ChickWeight.\"\n",
    "\n",
    "The variables in this data include:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(ChickWeight)\n",
    "str(ChickWeight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that there are 4 variables - weight (numerical DV), Chick (observation ID), Time (number of days since birth when measurement was obtained - repeated measures IV), and Diet (4 level factor of different types of feed, categorical IV).  If we look at some observations, we'll see that there are multiple repeated measures for the weight of the same chicken at different time points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(ChickWeight, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at just the first 10 rows, we see that all of these rows refer to the same chick (id = 1) at multiple time points.  If you recall in our one-way repeated measures example, we had wide data (the time points were all in the same row, different columns) but here our data is already in \"long\" format so we don't have to do that bit of data cleaning.\n",
    "\n",
    "We run this analysis with the same setup as all other analyses, with the time as it's own variable (in this case it's treated as categories and not numbers).  But there is one slight difference, the model is specified as:\n",
    "\n",
    "weight~ (Diet*Time) + Error(Chick/(Time)\n",
    "\n",
    "You'll notice the addition of an Error term. This error term is required because of the multiple observations arising from the same units (in this case Chicks).  This is the part of the residual variance (SS within) that is attributable to the individual unit (between unit variance) - such as genetics - some Chicks may just grow differently than others and it is not attributable to diet, but it is also not purely residual variance either.  This is the \"leftover\" that is explained by the Chick, and not by the predictors.\n",
    "\n",
    "Before we run the analysis, lets check our assumptions <a id=\"setup\"></a>\n",
    "### Repeated Measures - Assumptions Pre-Check\n",
    "\n",
    "First, I want to graphically display the means by time and diet to get a feel for the data.  We do this the same way as above, only this time one of our IVs is time and the other is a grouping variable (diet).  We'll make time our x variable and diet our color/line differentiator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## THIS OPTIONS IS ONLY FOR BINDER DISPLAY!!!\n",
    "options(repr.plot.width=5, repr.plot.height=4)\n",
    "\n",
    "#means by IVs (diet and time)\n",
    "ChickWeight %>% ggline(x = \"Time\", y = \"weight\", color = \"Diet\",\n",
    "       add = c(\"mean_se\", \"jitter\"),\n",
    "       size = 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even before we run our ANOVA we can see that we might have an interaction effect - after time 12 Diet 3 overtakes Diet 4 in highest mean weight, however the differences before time 12 between the Diets are so small that it may not be significant.\n",
    "\n",
    "Our other \"pre-flight\" checks:\n",
    "\n",
    "- DV is numeric (interval or ratio) - yes\n",
    "- No extreme outliers - nothing looks too bad on the graph, although I can see one chick in the green Diet that has extremely low weights after time 8 - did he stop growing?\n",
    "- Normality of residuals - can't check this until after\n",
    "- Independence of Observations (random selection, different samples) - We're going to account for the repeated nature of our measurements (multiple measurements in the same chicks) with our added Error() term.\n",
    "- Group sample sizes are approximately equal\n",
    "\n",
    "AND....\n",
    "\n",
    "Homogeneity of Variance - Let's check this right now with Levene's test.\n",
    "Recall:\n",
    "\n",
    "H0:  The variances in the groups are equal.\n",
    "HA:  The variances in the groups are not equal.\n",
    "\n",
    "In this test, we sort of want to fail to reject null, because it's easier if our variances are equal and we don't need to make the adjustment.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LeveneTest(DV ~ IV, data = your data frame)\n",
    "\n",
    "LeveneTest(weight ~ Diet, data = ChickWeight)\n",
    "LeveneTest(weight ~ as.factor(Time), data = ChickWeight)\n",
    "LeveneTest(weight ~ Chick, data = ChickWeight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"repeatanova\"></a>\n",
    "The variances are not equal, but we'll proceed with aov() anyway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ChickWeight$Time <- factor(ChickWeight$Time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 <- aov_car(weight~ (Diet*Time) + Error(Chick/Time), data = ChickWeight)\n",
    "knitr::kable(nice(a1))\n",
    "summary(a1$aov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows us our familiar output of our two main effects and our interaction effect, however it has adjusted the denominator of the F values for the within-Chick variance.  \n",
    "<a id=\"rsquarerep\"></a>\n",
    "We also get our effect size in this table as well - the \"ges\" is the generalized eta squared value.  It is not surprising that Time explains 77% of our variance, since that is the main thing that causes the chicks to grow - we assume that they will grow over time, it would be bad if they didn't.  What is really important here instead are the effects of Diet and the interaction of Diet over Time.  As we saw in the graph, the Diet itself appeared to have more of an effect on the weights at each time period than the overall trajectory of weights over time (the interaction).\n",
    "\n",
    "Compare this to an ANOVA without the Error() term.  This time we don't have two sets of Residuals, only one.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(aov(weight~ (Diet*as.factor(Time)) , data = ChickWeight))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"residnormrep\"></a>\n",
    "### Post-Hoc Assumptions Check - Normality of Residuals\n",
    "As always, we need to do our post hoc checks for normality of residuals. The basic \"aov\" object is nested inside the aov_car object, so we take that into account in our call to plot.  Also we have two sets of Residuals to inspect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## THIS OPTIONS IS ONLY FOR BINDER DISPLAY!!!\n",
    "options(repr.plot.width=5, repr.plot.height=4)\n",
    "qqnorm(proj(a1$aov)$Chick[,\"Residuals\"])\n",
    "qqline(proj(a1$aov)$Chick[,\"Residuals\"])\n",
    "\n",
    "qqnorm(proj(a1$aov)[[3]][,\"Residuals\"])\n",
    "qqline(proj(a1$aov)[[3]][,\"Residuals\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}