content/lab/07_category_fx/lab07_blank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../../shared/img/banner.svg\"></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 07 - Modeling Categorical Effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "sys.path.append(\"../../\")\n",
    "\n",
    "from shared.src import quiet\n",
    "from shared.src import seed\n",
    "from shared.src import style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "from client.api.notebook import Notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import scipy\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shared.src.utils.util as shared_util"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_context(\"notebook\", font_scale=1.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ok = Notebook(\"ok/config\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "1. Practice writing categorical-numerical models in pyMC3\n",
    "1. Become more familiar with running and reporting the results of classical statistical tests with `scipy.stats`.\n",
    "1. Recognize the benefits of the modeling flexibility of the sampling approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ❤️ React for Bayesian Monte Carlo, 😡 React for Null Hypothesis Significance Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You're working as a research scientist studying\n",
    "social cognition and social processing.\n",
    "You're using fMRI machines\n",
    "to search for areas of the brain that respond to emotionally and socially salient stimuli.\n",
    "\n",
    "Since your grant money comes from Facebook,\n",
    "you're trying to find areas of the brain\n",
    "that respond differentially to discovering that\n",
    "a social media post has received one of the six reactions:\n",
    "`like`, `love`, `haha`, `wow`, `sad`, and `angry`,\n",
    "pictured below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ![reactions](https://upload.wikimedia.org/wikipedia/en/thumb/4/49/Facebook_Reactions.png/309px-Facebook_Reactions.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From a literature review, you've identified a region that is likely to respond to these stimuli.\n",
    "\n",
    "You've collected the average activation values for the center of this region\n",
    "for a number of participants as they observed reactions to their social media posts.\n",
    "The data is stored in `fmri_data.pkl`. (This data is simulated, not real)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_df = pd.read_pickle(\"data/fmri_data.pkl\")\n",
    "\n",
    "fmri_df.sample(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_df.reaction.unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The area of interest takes up one \"pixel\" of the fMRI image.\n",
    "Since fMRI measures brain activity in 3D space,\n",
    "that is, in a volume,\n",
    "these are known as _voxel_ s,\n",
    "short for \"volumetric pixels\".\n",
    "The values you've measured are called _voxel activities_\n",
    "or _voxel activations_,\n",
    "shortned to `voxel_act` in the `fmri_df`.\n",
    "\n",
    "For each participant,\n",
    "you also measured the average activation of this voxel during the presentation\n",
    "of non-socially-relevant, non-emotional stimuli,\n",
    "and subtracted this value off.\n",
    "\n",
    "Therefore, you're looking to see whether the average activation across subjects (rows of `fmri_df`)\n",
    "is different from 0.\n",
    "You suspect that an ANOVA is the right model for this data:\n",
    "a categorical effects model with a Gaussian likelihood\n",
    "and equal variances in each group."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Begin by visualizing the data.\n",
    "\n",
    "Consider a box-and-whisker plot, or a collection of distplots.\n",
    "\n",
    "#### Q Does this look like data on which this model can be used? Explain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define an ANOVA model for this data in pyMC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define priors for the group means,\n",
    "making sure there is a mean for each group,\n",
    "and then use the group labels (`react_idx`)\n",
    "to set the parameters of a `Normal`ly-distributed variable\n",
    "for the likelihood.\n",
    "For that variable, set `observed` equal to the observed voxel activations.\n",
    "You'll also need a prior for the (pooled) standard deviation.\n",
    "You are free to choose your prior as you want.\n",
    "\n",
    "Refer back to the lecture material if you get stuck."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Explain your choices for the prior component of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using this model, draw samples and visualize the posterior with `pm.plot_posterior`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Which of the group means are different from 0? What does that mean for your original hypothesis that this brain area responds to emotionally-salient stimuli?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, run a traditional ANOVA on the data using `scipy.stats.f_oneway`.\n",
    "\n",
    "Save the values of $F$ and $p$ to `reaction_F` and `reaction_p`, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ok.grade(\"q1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Was the difference in group means statistically significant? Does this agree or disagree with your results from the pyMC model above?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beyond ANOVA: Modeling 👏 Emoji 👏 Use 👏"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not every kind of effect can be detected by an ANOVA,\n",
    "since not all data is normal and not all changes affect the mean.\n",
    "\n",
    "Say you're working in a lab studying [Human-Computer Interaction](https://en.wikipedia.org/wiki/Human%E2%80%93computer_interaction).\n",
    "\n",
    "Your thesis is on emoji use.\n",
    "What causes people to use emojis?\n",
    "How can we improve the emoji experience for users?\n",
    "\n",
    "As an early step in this research, you're collecting some demographic data on emoji use.\n",
    "You suspect that younger generations (those born after 1985, or with a birthdate `>1985`) use more emojis than do older generations (those born before 1985, or with a birthdate `<1985`).\n",
    "To test your hypothesis, you collect some data, `emoji_data.pkl`. (This data is simulated, not real)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emoji_df = pd.read_pickle(\"data/emoji_data.pkl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emoji_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `emoji_ct` variable counts the number of emojis in the text message,\n",
    "while the `generation_idx` and `birthdate` variables code the generation of the text message sender,\n",
    "the former as a number and the latter as a string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emoji_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the distribution of the `emoji_ct` for each generation.\n",
    "\n",
    "#### Q Use your visualization to predict whether an ANOVA will return a positive result. Explain your prediction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Run an ANOVA and report the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result doesn't accord with your understanding of the data, so you decide to dig deeper.\n",
    "\n",
    "#### Plot a histogram of the data from each group, using `bins` as defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.arange(0, max(emoji_df[\"emoji_ct\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data within each group is decidedly not normally distributed.\n",
    "\n",
    "#### Q What's the most salient non-normal aspect of the data? Hint: Look especially at the `<1985` histogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the negative ANOVA result, the two distributions look different.\n",
    "\n",
    "#### Q Describe, in your own words, the differences between the two distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears that there are two kinds of text messages:\n",
    "messages _with_ emojis, which then have a variable number of emojis,\n",
    "and messages _without_ emojis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that we are working with counts,\n",
    "the natural distribution for modeling this data\n",
    "is the `ZeroInflatedPoisson` distribution\n",
    "(see [the pyMC docs for a plot of the distribution](https://docs.pymc.io/api/distributions/discrete.html#pymc3.distributions.discrete.ZeroInflatedPoisson)).\n",
    "\n",
    "You can think of it as a random variable with two parts:\n",
    "its value is either equal to 0, or it is Poisson-distributed.\n",
    "You could implement it from scratch with `pm.math.switch`.\n",
    "It has two parameters:\n",
    "`theta`, which is the rate of the Poisson part,\n",
    "and `psi`, which is the chance the value is Poission-distributed.\n",
    "\n",
    "In our case,\n",
    "`psi` is the chance that the text message contains any emojis,\n",
    "while `theta` is the average number of emojis in the texts that contain emojis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Explain, in your own words, why a sampling or Monte Carlo approach, like the one we are taking in this class, allows us to work with zero-inflated Poisson random variables more easily than a traditional analytical approach, based on defined statistical tests like the $F$ test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, define a model for the emoji data that uses the Zero-Inflated Poisson.\n",
    "\n",
    "Write it by analogy to the ANOVA model:\n",
    "define priors for the model parameters,\n",
    "making sure each group gets a random variable\n",
    "for each of its prameters\n",
    "(notice that the `ZeroInflatedPoisson`\n",
    "variable has two parameters, not just one,\n",
    "and neither is shared, unlike the `sd` in an ANOVA).\n",
    "Then index into those with the group labels (`generation_idx`)\n",
    "to set the parameters of the `ZeroInflatedPoisson` variable.\n",
    "For that variable, set `observed` equal to the observed counts.\n",
    "\n",
    "Think carefully about the priors for the parameters of the Zero-Inflated Poisson variable.\n",
    "There are multiple right answers, but also lots of wrong answers.\n",
    "The variable `psi` must be between `0` and `1`\n",
    "and is continuous,\n",
    "while the variable `theta` just has to be above `0` (positive-only).\n",
    "You might have some other prior information or assumptions you want to include\n",
    "by choosing the initial distribution for these variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (2 pts) Using this model, draw samples and visualize a posterior over the parameters `psi` and `theta` with `pm.plot_posterior` or something similar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Q Report your findings. Do the groups differ? If so, how?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ok.score()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}