1 - repeated_measures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "efcfa853-a47e-46e7-952d-e6ca034e602e",
   "metadata": {},
   "source": [
    "# Working with Repeated Measures Data\n",
    "---\n",
    "\n",
    "## Learning Objectives\n",
    "- Learn how to use split-apply-combine to average repeated measures data appropriately\n",
    "- Learn how to properly visualize repeated measures data with Seaborn\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "974eda74-5a16-4874-8c30-a20d446fe820",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In the previous chapter on Seaborn plots, we worked with some simulated reaction time (RT) data from 5 individual participants. We showed various Seaborn plots that allowed us to examine the variance in the data sets. However, for the sake of focusing on the different types of Seaborn plots and what they show us, we glossed over one really important aspect of the data — its *nesting structure*, and in particular the fact that we have repeated measures. Repeated measures data are very common in cognitive psychology and neuroscience. In fact, while they are often treated as a \"special case\" in statistical and data science texts, in this context we focus on them because most of the data we will encounter later in this course include repeated measures — and **it's very easy to treat the data incorrectly and draw incorrect conclusions if you're not used to working with repeated measures data**.\n",
    "\n",
    "Recall that in the *Thinking About Data for Plotting* lesson, we learned that variance within individual units of measurement (such as individual participants) is not the same as the variance between participants. Furthermore, in a typical experiment our interest is in examining the average effects across individuals, and the range of differences between individuals. When we have nested data, such as many individual trials for each participant, this is because we recognize that any individual measurement may not be accurate, and so by averaging over many measurements we are likely to get a more accurate estimate of the true measurement for each individual. However, if we average the trial-to-trial variance of each individual (which is typically low, within the individual) across individuals, we lose the true inter-individual differences. \n",
    "\n",
    "For this reason, when plotting nested data it is usually appropriate to first average within each individual (for each condition), and then compute and visualize measures of the true inter-individual variability. Seaborn does not do this by default, and so the plots we generated in the previous lesson were actually inappropriate for assessing whether there was a true experimental effect. In this lesson we'll see how to use pandas `.groupby()` to use Seaborn appropriately with nested data. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb79b9d7-b000-4778-9936-69c54e0a73c2",
   "metadata": {},
   "source": [
    "### About the data\n",
    "First, let's load in a new data set (`flanker_rt_data.csv`). This real data from a flanker experiment (similar to that described in the section on [spreadsheets](../2/spreadsheets.html)).  In this experiment, participants had to press either the left or right arrow key, to indicate whether an arrow shown on the screen is pointing left or right, respectively. However, the catch is that the centre arrow is \"flanked\" by two other arrows on each side; these can be pointing the same way as the target arrow (congruent):\n",
    "\n",
    "![flanker_congruent](images/flanker_congruent@0.75x.png)\n",
    "\n",
    "or in the opposite direction (incongruent)\n",
    "\n",
    "![flanker_incongruent](images/flanker_incongruent@0.75x.png)\n",
    "\n",
    "We expect that the incongruent condition would be associated with slower RTs, because the flanking arrows create some visual confusion and response competition (cognitive interference; this process is not unlike that underlying the famous [Stroop](https://en.wikipedia.org/wiki/Stroop_effect) experiment, which you are likely familiar with). In order to confirm this hypothesis, we would want to get an estimate of the RT for each condition, from a representative sample of human participants. We need a number of participants because we know that there is variability in the average RT from person to person (**between-subject** variability). Here we have 27 participants.\n",
    "\n",
    "To estimate the average RT for each condition for an individual, we would want to present many trials of each condition, in random order. We do this because there is always some trial-to-trial variability in measuring human RTs (**within-subject** variability). In this study each participant completed 40 trials in each condition. Thus we have repeated measures of RT for each participant, in each experimental condition. \n",
    "\n",
    "This experiment also had a third condition, *neutral*. On neutral trials, only the centre arrow was presented, without flankers. This condition was included to control for the fact that 5 arrows all pointing in the same direction might facilitate (speed up) responses relative to a single arrow — there might be additive, beneficial effects of congruent flankers. Comparing both the congruent and incongruent conditions to the neutral condition allows us to estimate the relative cost or benefit of flankers that are incongruent or congruent with the direction of the centre arrow. To ensure that there were as many neutral (single-arrow) as flanker (five-arrow) trials, 80 neutral trials were included per participant. \n",
    "\n",
    "So we have data from 27 participants, and approximately 160 trials per participant in total.  The nesting structure of this data set is thus trials (n=160) nested within conditions (n=3) nested within participants (n=27). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "596938a1-bbb1-4e4e-b992-1e9521465f88",
   "metadata": {},
   "source": [
    "## Load the data\n",
    "\n",
    "First we'll import the packages we'll need, read in the data, and view a sample of 10 rows of the DataFrame:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "332ca78c-05e2-4544-9555-50365b324b05",
   "metadata": {},
   "source": [
    "~~~python\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "df = pd.read_csv('data/flanker_rt_data.csv')\n",
    "\n",
    "df.sample(10)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c1fd932-ae73-4021-89ac-299cf29b79fd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "41fdcd00-8be7-43d2-aa55-c6ef513bb32b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Examine the data without considering repeated measures\n",
    "\n",
    "First we'll get some descriptive statistics for RT in each condition. We do this using pandas `.groupby()` method to group by flanker condition, and then we select the `rt` column for further processing by putting it in square brackets (which is pandas' way of selecting columns). Finally we chain all this with the `.describe()` method:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11a6e56-50dd-4e85-98d6-7fd2d063914a",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "df.groupby('flankers')['rt'].describe()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa35f398-ec2e-4492-9bb3-c61e7cb77573",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ebf678ca-246e-4835-b601-d64c2d03cff5",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "Note the `count` (total number of data points), standard deviation (`std`) as well as the `min` and `max` values for each condition.\n",
    "\n",
    "Let's look at a box plot of values in each condition:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e709200-b36a-4ab9-819e-abb47e3f41e3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='box',\n",
    "           data=df,\n",
    "           x='flankers', y='rt')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76214dd5-0f0d-43c0-89fe-b41cfe6a2b68",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c3d2e43d-ccf5-485b-a345-18def5cf475b",
   "metadata": {},
   "source": [
    "The data show the skew that is expected of RT data — there are fundamental human performance limits on how fast a person can response (short RTs), but the upper limit of RTs is theoretically unlimited. In the present experiment, participants actually only had a 1 s window in which to respond, so the RTs are all less than 1. Nevertheless, the data are skewed because there is a wider range of RTs above the median than below. \n",
    "\n",
    "Another thing we can see about the distribution of values is that there are many outliers (individual points in the box plot) at the long end of the RT range, but few at the short end. This is, again, a consequence of the skewed nature of the data. However, notably some of the outliers on the short end of the distributions are suspicious: human RTs in very simple tasks (like pressing a button any time a light comes on) are rarely shorter than 200 ms, yet we have a few RT values below that. So our box plots suggest that we may have some anomalous data, perhaps from trials on which someone pressed a button accidentally prior to processing the actual direction of the arrow.\n",
    "\n",
    "Finally, the box plots suggest that RTs are slowest in the incongruent condition and fastest in the congruent condition — consistent with our predictions based on past studies of the flanker effect. RTs for the neutral condition are in between the other two conditions. This suggests that perhaps both facilitation effects in the congruent condition, and interference effects in the incongruent condition, are present in the data. However, it's hard to know how large or believable those differences are, because the distributions overlap quite a bit between the three conditions.\n",
    "\n",
    "### Stingray plot\n",
    "Let's look at a couple of other ways of viewing the data that we learned about in the Seaborn lesson. First, we'll use a violin (stingray) plot to see both the KDE and box plot:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02d3ad90-f871-4c97-b2ab-4eb698a3af7b",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='violin',\n",
    "           data=df,\n",
    "           x='flankers', y='rt')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e6a3a3e-033d-42a9-a930-1550b51ff7c1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5ef2a1d8-62fb-4438-a178-eedd7ee74843",
   "metadata": {},
   "source": [
    "The skewed nature of the data are again evident, but we can see that not only are the median values for the three conditions different, but the bulk of values in each distribution of RTs also appear to be rather different. That is, the \"humps\" in the middle of each distribution are relatively non-overlapping, although neutral and incongruent are fairly similar.\n",
    "\n",
    "### Strip plot\n",
    "Next let's look at a strip plot of the data, which plots each data point separately:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60adbb71-2dde-458e-bbe4-782112874a61",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='strip',\n",
    "           data=df,\n",
    "           x='flankers', y='rt')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e8b848-2401-4e2b-a029-919aee739e0f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "04db5de3-0d57-4ab6-b5a9-db8767aef852",
   "metadata": {},
   "source": [
    "This is really less informative than the above plots, because there are so many data points! Mostly we just see big patches of the plots for each condition with lots of individual values. \n",
    "\n",
    "However, the strip plots do show us something interesting that wasn't evident in the previous plots: the very fast outlier data points are quite anomalous, even compared to the very slow outliers. That is, there are few \"fast\" data points, and they are generally quite far from the rest of the distribution of RTs. In contrast, the slower data points are more plentiful, and seem to form a rather continuous set of values, albeit more sparse than the points closer to the centres of the distributions. This suggests that the slower data points might be quite believable as real data — mosly people make faster resposnes, but sometimes they're a bit slower. In contrast, the strip plot gives us greater confidence that the very fast RTs should be eliminated from analysis, because they are not just slightly-faster reponses than he next-fastest responses.\n",
    "\n",
    "### Point plot\n",
    "Finally, let's look at a point plot. Recall that point plots show the mean for each condition, along with 95% confidence intervals (CIs) as the measure of variance. CIs are much narrower than the whiskers in box plots, because they are meant to represent our confidence in the estimate of the mean, rather than showing the overall distribution of values. In general, we can interpret conditions with non-overlapping 95% CIs as being \"believably different\" from each other, and likely to be statistically significantly different as well."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d1b7657-75fb-4609-90b8-5e7e0b1373a4",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='point', join=False,\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='flankers')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8b80039-bf6e-4dab-8b99-fb68481ff845",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0a1739fc-2f6c-4a3e-8bc1-4b136f3ae509",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.set_context('talk')\n",
    "sns.catplot(kind='box',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', col='participant', col_wrap=6,\n",
    "           )\n",
    "plt.show()\n",
    "sns.set_context('paper')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8296754-5e9a-4151-8b63-ab584e77af87",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "647be30d-be3e-4cc9-a783-a7232090f112",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Examine individual participant data\n",
    "Now, let's use split-apply-combine to average the RTs within each individual participant (and within each condition within individuals) prior to the group-level descriptive statistics. We can do this in one line of code by chaining, combining the following steps:\n",
    "- group by both participant ID and condition (split). Note that the argument to `.groupby()` is a list containing the two columns we want to group by\n",
    "- compute the mean for each condition within each participant (apply)\n",
    "- group by condition only (combining across participants)\n",
    "- use `.describe()` to apply another computation (average across participants) and combine the output into a new table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0e1434-bf44-4862-afae-3e592229f67f",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='strip',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='participant')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b30d6dd2-b3b5-4f06-899a-e5e2abfe9ac3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "168612e2-f905-42fe-8505-5921c4c9f21b",
   "metadata": {},
   "source": [
    "This is a lot of numbers to look at! Indeed, so may that pandas doesn't show us all the rows. From what is shown, we can see that there is certainly variance in mean RTs for each individual, and in variance (as measured by standard deviation), but it's hard to get a real sense of how much variability there is, or whether there are any participants whose data may be anomalous relative to most others.\n",
    "\n",
    "### Visualizing nested data with Seaborn\n",
    "\n",
    "Seaborn makes it easy to separate data by participants in visualization. For example, we can generate box plots as we did above, but put participant on the *x* axis, and use the `hue` kwarg to colour-code flanker condition. Because there are lots of individual participants, we add the `aspect=3` kwarg to tell Seaborn to make the plot 3x as wide as it is tall."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18774228-dbc5-4a39-830d-e76f9c5a10ff",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='swarm',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='participant')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e80aa32-e974-4f2b-a6a3-49fcbe6165d4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "445f1f6e-ce91-428f-83ff-b9e65a7d00d3",
   "metadata": {},
   "source": [
    "Another way to visualize so many different box plots is with Seaborn's `col` and `col_wrap` kwargs. On its own, `col=` will plot whatever variable you give it in separate columns — in our case, participants (rather than as different locations on the *x* axis, as shown above). Adding the `col_wrap` kwarg is necessary when the number of categories passed to `col` is large, as in our 27 participants. This kwag tells Seaborn how many columns to put on one row before wrapping around to a new row. \n",
    "\n",
    "One last trick we will use is to set Seaborn's context to `talk`, to make the font size easier to read when we have so many plots. We'll set it back to `paper` after showing the plot, so our future plots aren't affected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e106d5cc-1ff5-4e72-8ffb-96e524b2aef5",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "~~~python\n",
    "\n",
    "sns.set_context('talk')\n",
    "sns.catplot(kind='box',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', col='participant', col_wrap=6,\n",
    "           )\n",
    "plt.show()\n",
    "sns.set_context('paper')\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "050b57fc-af60-49f1-bdf1-3cb6c7263213",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2cfb3746-c904-43f6-90d7-3f18f2d129fa",
   "metadata": {},
   "source": [
    "This plot shows us that, indeed, there is variability among participants in both the median RTs and the variance (as indicated by the length of the whiskers, and the number and range of outlier points). Encouragingly, however, each individual seems to show the expected effect, of longer RTs for incongruent than congruent flakers. It also shows us that those anomalously fast RTs we saw above are contributed only by a couple of participants, again leading us to think that these are not believable data points.\n",
    "\n",
    "Next we'll generate a strip plot, in which we colour-code the points according to which participant they are from:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd4311d2-ef66-4b76-8c1b-ed38d596f5d2",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='strip',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='participant')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6413300e-61b7-47e4-9572-e3ff576947f6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "39a77877-9651-445a-9c86-13d23aa5f982",
   "metadata": {},
   "source": [
    "This is not a lot more useful than the strip plot above, because Seaborn draws dots for each participant over those of previous participants. So we see lots of pink in the plot above, because that is the colour used for the last participant plotted. This is addressed by using a swam plot (which will take a long time to execute if you try to run this, due to the large number of data points it's trying to plot):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bbe881d-fc73-4560-856e-a51f8acae8a1",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='swarm',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='participant')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48ff4515-0790-4e6c-b340-db936df4d045",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3965214f-7d7e-4b8e-8957-4f2b6a1c45d4",
   "metadata": {},
   "source": [
    "However, in this case there are so many data points that only a relatively small percentage of them are actually drawn (hence the warning messages above). It's hard to get a sense of the distributions of values for individual participants. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47e4c3b1-439e-443b-a372-5866297d8b1e",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Appropriately dealing with repeated measures: Aggregating data by participants\n",
    "\n",
    "Because there are so many data points, and variation in both the mean and range of values between participants, it is necessary to first average (**aggregate**) within each participant and condition, prior to computing descriptive statistics or visualizing. This ensures that the variance we are examining is the **between-subjects variance**, unconfounded with **within-subjects variance**. This is primarily important in generating plots that show confidence intervals (CIs), such as point plots and bar plots, for reasons explained below.\n",
    "\n",
    "Seaborn has a way of doing this \"on the fly\" when you generate plots, using the `units=` kwarg and passing the name of the grouping variable for repeated measures (in this case, `participant`):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c0e3b32-82f0-45bb-b46b-864250bbec40",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='point', join=False,\n",
    "           data=df, units='participant',\n",
    "           x='flankers', y='rt', hue='flankers')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9de7dab0-d0de-4f4e-a04f-22006ccc390a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9ef23619-3572-49ce-b300-ec537518570d",
   "metadata": {
    "tags": []
   },
   "source": [
    "**This is the most important figure in this lesson**. Compared to the point plot with 95% CIs that we generated above for the aggregated data, the CIs in this plot are much larger than those we saw earlier, in the point plot for the un-aggregated data. In that previous plot, the width of the CIs was only about 0.005 s (5 ms). Here, with the CIs computed properly (as reflecting the variance across individuals), the width of the CIs is approximately 0.02 s (20 ms). While still small relative to the differences in the means, the CIs are nonetheless roughly 4x wider in the aggregated data. \n",
    "\n",
    "### Why do CIs get larger when we properly aggregate repeated measures data?\n",
    "\n",
    "We've already seen that the variance in the data actually went down, as quantified by the min-max range of values, and the standard deviation. So it may seem strange that the CIs do the opposite, and get bigger. To understand this, consider how the 95% CIs are computed:\n",
    "\n",
    "$CI = 1.96 * \\frac{std}{\\sqrt n}$ \n",
    "\n",
    "CIs are computed by dividing the standard deviation ($std$) by the square root of $n$ — the number of data points used to compute the mean and standard deviation. In the un-aggregated data, this was approximately 1000 – 2000 data points for each condition. In contrast, for the aggregated data, $n$ = 27. Thus the CIs are much wider, because dividing by a larger number (1000 - 2000, rather than 27) results in a smaller number (even factoring in the differences in standard deviation, and the square root calculation). The values are multiplied by 1.96 because this is the *z* value corresponding to *p* = .05 in a normal distribution.\n",
    "\n",
    "Notably, the conclusions that we would draw from this final plot would be quite different from if we had, in error, used the first point plot generated from the un-aggregated data. The first plot made it look like RTs in the three conditions were very different from each other, with little question as to the differences between them (due to the small and non-overlapping CIs). In contrast, when the data are properly aggregated as in the above plot, the CIs for neutral overlap a lot with both of the flanker conditions, which makes us less confident in their differences. \n",
    "\n",
    "\n",
    "Recall that the neutral condition was intended to help assess whether the congruent flankers showed facilitation relative to neutral, and whether incongruent flankers showed interference relative to neutral. The fact that the neutral condition overlaps with both others, but congruent and incongruent appear significantly different, suggests that the flanker effect (congruent faster than incongruent) may actually be due to a *combination* of facilitation and interference, rather than one or the other. \n",
    "\n",
    "The interpretation of the current data set is less important here than understanding the importance of recognizing and appropriately dealing with repeated-measures data. It is critical to understand the effects of repeated measures and nesting on estimate of variance — and of the differences between measures of variance in the data (as reflected by things like standard deviation, or box plots), and estimates of our confidence in the accuracy of the mean (CIs). CIs can be used for statistical inference (i.e., is there a believable difference?), whereas the variance reflected in box plots (and histograms, KDEs, etc.) is more informative at the early stages of exploring and understanding your data, and diagnosing potential issues in the data that might need to be addressed prior to inference about the true patters in the data. Although this initial example used relatively \"clean\" data — and we didn't worry about doing anything to address the skew of the distributions or the outliers — in future lessons we will see cases where this first step can be critically important. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba2606b8-fc34-441e-a8e3-2fd202273792",
   "metadata": {},
   "source": [
    "## Aggregating using Split-Apply-Combine\n",
    "\n",
    "While it's nice that Seaborn can aggregate repeated-measures data for us, that only works for generating Seaborn plots. It's also good to know how to create a DataFrame of aggregated data, which you can then use to generate tables of descriptive statistics, and eventually perform statistical tests on.\n",
    "\n",
    "To do this, we will use *two* split-apply-combine operations in sequence. In one chained command, we will:\n",
    "- split the data by participant and condition: `.groupby(['rt', 'flankers'])`\n",
    "- select the `rt` column: `['rt']`\n",
    "- apply a calculation to get mean RT for each participant/condition: `.mean()`\n",
    "- split these results by condition only, to combine across participants: `.groupby('flankers')`\n",
    "- apply descriptive statistics and combine the results in a table: `.describe()`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f552407-4b9a-47be-98b1-ef198053a1ec",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "df.groupby(['participant', 'flankers'])['rt'].mean().groupby('flankers').describe()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b35412c4-b0dc-47cf-a10b-fafa993ede1f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0809a111-ab72-451b-bc95-6d60fd21d887",
   "metadata": {},
   "source": [
    "Compare the `count`, `std`, `min`, and `max` values here to that in our first descriptive statistics table above, which we'll reproduce here for ease of comparison:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44f4b167-1869-4e5f-939c-bc4dcc3191e3",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "df.groupby('flankers')['rt'].describe()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6738b9b-08b8-4362-8916-0c05328cfd78",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "90d4fd05-62c9-489c-92ba-612dce817659",
   "metadata": {},
   "source": [
    "You can see that the number of data points goes down from the total number of trials across all participants (~4200) to 27 per condition (the number of participants). The mean RTs do not actually change at all. This is expected because averaging is a linear operation, meaning that the average across all trials is the same as averaging first within participants then across trials. \n",
    "\n",
    "However, the range of values decreases when we look at the `min` and `max` values. This is a normal consequence of averaging — the most extreme values found in a set of individual data points will be reduced when we average by some grouping variables (like participant).  \n",
    "\n",
    "The variance (as reflected by the standard deviation, `std`) gets smaller as well, for a similar reason: the variance is now among the average RTs for each participant in each condition, which have a smaller range than the raw data points. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "297806bf-0a31-4b94-a6c2-112f95bdacd8",
   "metadata": {},
   "source": [
    "### Plot data aggregated across trials\n",
    "\n",
    "Now let's re-plot the data using the same plots we did above, but now for the data averaged within participants. \n",
    "\n",
    "First let's save the aggregated data in a new DataFrame called `df_avg`. We chain the `.reset_index()` method to the end of this, because `.groupby()` sets the indexes of the DataFrame to be the grouping variables (in this case, `rt` and `flankers()`). This is ideal fo rlooking at tables of the data (e.g., in the above, RT and condition are in bold), but index columns are not visible to Seaborn as data columns for plotting. So we use `.reset_index()` to convert `rt` and `flankers` back from indexes to regular DataFrame columns."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b4653eb-f3cd-4061-9389-7feedfbda1aa",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "df_avg = pd.DataFrame(df.groupby(['participant', 'flankers'])['rt'].mean()).reset_index()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eba0b0d7-1b88-4cf9-82c2-3529e92b76b9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ddf961c2-9bc6-4707-89d0-28ad6eb8a381",
   "metadata": {},
   "source": [
    "Now let's generate a box plot of these data:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b07cf4a-84ce-4070-adcb-8af925999929",
   "metadata": {},
   "source": [
    "~~~python\n",
    "sns.catplot(kind='box',\n",
    "           data=df_avg,\n",
    "           x='flankers', y='rt')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcfba48d-96a4-4f36-ab49-8476bc330470",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f895f14d-3849-4051-9f71-0df2d087beec",
   "metadata": {},
   "source": [
    "Another, somewhat confusing, thing is that the order of the conditions on the *x* axis is reversed: whereas before the order was congruent- neutral - incongruent, now neutral is last. This is a consequence of the `.reset_index()` method, which led to the conditions being re-ordered alphabetically. This highlights that its always good to pay close attention to your plot labels!\n",
    "\n",
    "Seaborn actually has an easy way to specify the order of values on a categorical axis, if you want to: using the `order=` kwarg, and passing a list of conditions in the order you want them to appear:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42411117-8032-402a-9c07-7590dcfa6a36",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='box',\n",
    "           data=df_avg,\n",
    "           x='flankers', y='rt',\n",
    "           order=['congruent', 'neutral', 'incongruent'])\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "877899e4-13d4-4463-bb4b-9bccef889555",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "451f5aa1-b4df-4b68-8f5d-52c47df12485",
   "metadata": {},
   "source": [
    "If you compare this to the box plot based on every data point that we generated earlier, you'll see that the range covered by the whiskers and the box are both smaller. This is consistent with the reduced range in the data that we observed in the tables of descriptive statistics. Likewise, the boxes (which represent the IQR) are now non-overlapping, whereas before they overlapped. \n",
    "\n",
    "### Strip plot\n",
    "A strip plot of the data now much easier to make sense of, because we only have one point per participant in each condition:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1040708-fc58-4c67-89ab-c274abd28ad5",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='strip',\n",
    "           data=df_avg,\n",
    "           x='flankers', y='rt', hue='participant')\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b92f7bff-1c08-4f3a-a6e8-c0a5fb47198c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fe7ffafb-d8cc-4b81-8f36-5e1793687186",
   "metadata": {},
   "source": [
    "### Line plot\n",
    "This makes it easier to get a sense of the variability between individuals. However, It's a bit hard to see how each individual's RT changes between the two experimental conditions. To facilitate this, let's use a point plot, which by default draws lines to join data from each individual:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356bcaab-f1dd-43d7-aca9-492739c9c1ee",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='point',\n",
    "           data=df_avg,\n",
    "           x='flankers', y='rt', hue='participant',\n",
    "           order=['congruent', 'neutral', 'incongruent'])\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3812f947-94b7-4614-bd43-efe0cbe10977",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4be1f5cd-65a1-48bb-a525-10fd7cdf9793",
   "metadata": {},
   "source": [
    "You may recall that point plots typically show means and 95% CIs, but the plot above lacks CIs. This is because in the `df_avg` DataFrame, there is only one value for each participant in each condition (because we averaged across trials). If we wanted to see CIs for each individual, we could use the un-aggregated data (`df`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cec359cc-3c42-4f4a-aaf2-249be3d5b846",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='point',\n",
    "           data=df,\n",
    "           x='flankers', y='rt', hue='participant',\n",
    "           order=['congruent', 'neutral', 'incongruent'])\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4b5681f-cc51-4c61-ba2c-d8882a81d2d4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b45550d6-e4a2-43f5-a18c-9857fc69a465",
   "metadata": {},
   "source": [
    "## Point plots aggregated across participants\n",
    "\n",
    "The CIs in the above plot don't, however, help us to make a simple inference as to whether, across all participants, there is a significant RT difference between congruent and incongruent. The data are certainly suggestive of this, because all participants show slower RTs for incongruent. But this is where a point plot with 95% CIs, across participants, is useful:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b943cc61-e9bf-4ef4-b224-a8d99d312abd",
   "metadata": {},
   "source": [
    "~~~python\n",
    "\n",
    "sns.catplot(kind='point', join=False,\n",
    "           data=df_avg,\n",
    "           x='flankers', y='rt', hue='flankers',\n",
    "           order=['congruent', 'neutral', 'incongruent'])\n",
    "plt.show()\n",
    "~~~"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f82293e-1dc1-4a7c-b3e4-e2c72754cd67",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e8fc3cf3-d510-481e-b79d-09db016b4836",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Summary\n",
    "- Repeated measures occur when we take repeated samples of the same measurement. This is common in cognitive psychology and neuroscience data, because individual measurements typically include noise. Averaging across repeated measurements yields a more robust estimate of the true value for that individual\n",
    "- When we have repeated measures data, we need to be concerned about the difference between within-subjects variance, and between-subjects variance\n",
    "- Within-subjects variance is the variation within an individual, across the repeated measurements we take\n",
    "- Between-subjects variance is the variation between individuals, in the *average* response across the repeated measures taken from that individual\n",
    "- If we simply average across all trials and participants, our measures of variance will be inaccurate, because we are conflating within- and between-subjects variance\n",
    "- Conflating within-and between-subjects variance will lead to incorrect conclusions, particularly because confidence intervals and statistical tests consider the number of samples in their computations. When assessing effects across a group of participants, the number of participants is the appropriate number of samples to use. The much larger number of samples that comes from summing the number of individual trials across all participants will artificially inflate estimates of the significance of experimental effects\n",
    "- The proper way to work with repeated measures data is to first average across all trials for each experimental condition, within an individual. We call this *aggregation*\n",
    "- After we have properly aggregated the data within each individual, then we can average across participants and accurately estimate the between-subjects variance\n",
    "- Between-subject variance is typically what is important in assessing the significance of experimental manipulations\n",
    "- Proper aggregation of repeated measures data is easy using the pandas `.groupby()` method"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  },
  "nteract": {
   "version": "0.28.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}