{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multilevel Modeling and Time-Series Forecasting\n",
    "\n",
    "In this lab you will fit multilevel (mixed-effects) models to nested behavioral data and compare them to a single-level regression that ignores the nesting structure. The dataset contains session-level response rates from eight participants, each observed across twenty sessions under varying reinforcement rates.\n",
    "\n",
    "**Learning objectives:**\n",
    "- Understand why ignoring nesting leads to biased inference\n",
    "- Fit and interpret mixed-effects models using `statsmodels`\n",
    "- Compare models using AIC and BIC\n",
    "- Visualize participant-level variation in intercepts and slopes\n",
    "\n",
    "**Required packages:** `pandas`, `numpy`, `matplotlib`, `statsmodels`, `scipy`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Run the cell below to import the libraries you will need."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Import libraries\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Load and Explore the Data\n",
    "\n",
    "Load `nested_behavior_data.csv` into a pandas DataFrame. The dataset has four columns:\n",
    "\n",
    "| Column | Description |\n",
    "|---|---|\n",
    "| `participant_id` | Participant identifier (1-8) |\n",
    "| `session` | Session number (1-20) |\n",
    "| `reinforcement_rate` | Reinforcements delivered per minute |\n",
    "| `response_rate` | Responses per minute |\n",
    "\n",
    "After loading:\n",
    "1. Print the first few rows and the shape of the DataFrame.\n",
    "2. Compute summary statistics grouped by participant.\n",
    "3. Create a scatter plot of `reinforcement_rate` vs. `response_rate`, coloring points by participant."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load the data\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Summary statistics by participant\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Scatter plot colored by participant\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Fit a Single-Level OLS Regression\n",
    "\n",
    "Fit an ordinary least squares (OLS) regression predicting `response_rate` from `reinforcement_rate`, pooling all participants together. This model ignores the nesting of observations within participants.\n",
    "\n",
    "1. Fit the model using `smf.ols()`.\n",
    "2. Print the summary.\n",
    "3. Plot the residuals vs. fitted values. Do you see any patterns that suggest the model is misspecified?\n",
    "4. Color the residual plot by participant. What do you notice?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Fit OLS model\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Residual plot colored by participant\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Why might the standard errors from the OLS model be misleading when the data are nested? Write your answer below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Fit a Random-Intercept Model\n",
    "\n",
    "Now fit a multilevel model that allows each participant to have their own intercept, while the slope of `reinforcement_rate` is fixed across participants. Use `smf.mixedlm()` from statsmodels.\n",
    "\n",
    "```python\n",
    "model = smf.mixedlm(\"response_rate ~ reinforcement_rate\",\n",
    "                     data=df,\n",
    "                     groups=df[\"participant_id\"])\n",
    "result = model.fit()\n",
    "```\n",
    "\n",
    "1. Fit the model and print the summary.\n",
    "2. What is the estimated variance of the random intercept?\n",
    "3. How does the fixed-effect slope compare to the OLS slope from Task 2?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Fit random-intercept model (use fit(reml=False))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Fit a Random-Intercept-and-Slope Model\n",
    "\n",
    "Extend the model so that each participant has both a unique intercept and a unique slope for `reinforcement_rate`. In statsmodels, you specify the random slope using the `re_formula` argument:\n",
    "\n",
    "```python\n",
    "model_rs = smf.mixedlm(\"response_rate ~ reinforcement_rate\",\n",
    "                        data=df,\n",
    "                        groups=df[\"participant_id\"],\n",
    "                        re_formula=\"~reinforcement_rate\")\n",
    "result_rs = model_rs.fit()\n",
    "```\n",
    "\n",
    "1. Fit the model and print the summary.\n",
    "2. Extract the random effects for each participant using `result_rs.random_effects`.\n",
    "3. For each participant, compute the participant-specific intercept and slope (fixed effect + random effect)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Fit random-intercept-and-slope model (use fit(reml=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Extract and display participant-specific intercepts and slopes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5: Compare Models Using AIC and BIC\n",
    "\n",
    "Compare the three models (OLS, random-intercept, random-intercept-and-slope) using information criteria. Lower values indicate a better trade-off between fit and complexity.\n",
    "\n",
    "**Note:** statsmodels does not populate `.aic`/`.bic` on `MixedLMResults` (they return `nan`), so we compute them from the log-likelihood. All models are fit by maximum likelihood here so the values are comparable.\n",
    "\n",
    "Create a table summarizing the AIC and BIC for each model. Which model is preferred?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Model comparison table (define info_criteria for the mixed models)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** What do the AIC/BIC values tell you about the importance of accounting for individual differences in this dataset? Write your answer below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 6: Interpret Fixed and Random Effects\n",
    "\n",
    "Using the random-intercept-and-slope model:\n",
    "\n",
    "1. State in plain language what the fixed-effect intercept and slope represent.\n",
    "2. State what the random-effect variances represent.\n",
    "3. Which participants have the steepest and shallowest slopes? What might this mean behaviorally (e.g., sensitivity to reinforcement rate)?\n",
    "4. Is there a correlation between participant intercepts and slopes? If so, what does it mean?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Analyze the relationship between participant intercepts and slopes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your interpretation here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 7: Visualize Participant-Level Regression Lines\n",
    "\n",
    "Create a figure with:\n",
    "1. A scatter plot of the raw data, colored by participant.\n",
    "2. The OLS regression line (single line for all data) in black.\n",
    "3. The participant-specific regression lines from the random-intercept-and-slope model, each in a different color.\n",
    "\n",
    "This visualization should make it visually clear why a single regression line is inadequate."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Visualization of participant-level regression lines\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 8: Visualize Random Effects\n",
    "\n",
    "Create two plots:\n",
    "1. A caterpillar plot (forest plot) showing the random intercepts for each participant with confidence intervals.\n",
    "2. A caterpillar plot showing the random slopes for each participant.\n",
    "\n",
    "**Hint:** You can extract the random effects from `result_rs.random_effects` and approximate confidence intervals using the conditional standard errors."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Caterpillar plots of random effects\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 9 (Optional): Time-Series Decomposition\n",
    "\n",
    "Select one participant and treat their session-by-session `response_rate` as a short time series.\n",
    "\n",
    "1. Plot the raw time series.\n",
    "2. Use `statsmodels.tsa.seasonal.seasonal_decompose` (with `model='additive'` and an appropriate `period`) to decompose the series into trend, seasonal, and residual components.\n",
    "3. Since our data may not have a true seasonal component over 20 sessions, focus on the trend and residual. What does the trend tell you about how this participant's behavior changed over sessions?\n",
    "\n",
    "**Note:** With only 20 observations, this decomposition is illustrative rather than definitive. The goal is to introduce the concept of decomposing temporal data."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Time-series decomposition for one participant\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflection\n",
    "\n",
    "In a few sentences, address the following:\n",
    "\n",
    "1. Why is multilevel modeling particularly important for behavioral research where data are collected across multiple participants and sessions?\n",
    "2. How would you explain the difference between fixed effects and random effects to a colleague who has only used standard regression?\n",
    "3. What are the limitations of the models you fit today (e.g., assumptions about the error structure, linearity)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your reflection here.*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}