Week10_TUT.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "7fbaf1e3",
      "metadata": {
        "id": "7fbaf1e3"
      },
      "source": [
        "# STA130 TUT 07 (Nov08)<br><br> 👪 📈 <u>Multiple Linear Regression</u>"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cdbec649",
      "metadata": {
        "id": "cdbec649"
      },
      "source": [
        "## ♻️ 📚 Review  / Questions [15 minutes]\n",
        "\n",
        "### 1. <u>Follow up questions and clarifications regarding concepts associated with **Simple Linear Regression**</u>\n",
        "    \n",
        "> First introduced on Oct21 Oct26, and conlcluded on Nov04 and Nov07... The theoretical \"normal distribution\" regression model, model fitting, hypothesis testing within this framework, and indicator variables and (two-sample) group comparision...\n",
        ">\n",
        "> - *You have GOT TO KNOW THIS ac rf bc*\n",
        "> - **Multiple linear regression** extends the **linear form** examined so far to multiple **predictor variables**; similarly extends the idea of **binary indicator variables** to **categorical variables**  with more than two levels; additionally provides both similar and extending **hypothesis testing** capabilities with the introduction of so-called **interaction variables** (which can be used, e.g., to examine the evidence that a linear association present between two variables differs across different groups of subpopulation groups in a dataset).\n",
        "\n",
        "<details class=\"details-example\"><summary><u><span style=\"color:blue\">The theoretical \"normal distribution\" regression model</span></u></summary>\n",
        "\n",
        "### The theoretical \"normal distribution\" regression model\n",
        "    \n",
        "$$\\Large Y_i = \\beta_0 + \\beta_1 x_i + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "- **Outcome** $Y_i$ is a **continuous numeric variable**\n",
        "    - **Outcome** $Y_i$ can also be called a **response**, **dependent**, or **endogenous variable** in some domains and contexts\n",
        "\n",
        "    \n",
        "- **Predictor variable** $x_i$ is a **numeric variable**\n",
        "    - Fow now we'll consider $x_i$ to be a **continuous** numeric variable, but this is not necessary, and we will consider other versions of $x_i$ later\n",
        "    - **Predictor variable** $x_i$ can also be called an **explanatory**, **independent**, or **exogenous variable**, or a **covariate** or **feature** (the last two of which are the preferred terms in the *statistics* and *machine learning* disciplines, respectively)\n",
        "\n",
        "    \n",
        "- **Intercept** $\\beta_0$ and **slope** $\\beta_1$ are the two primary **parameters** of a **Simple Linear Regression** model\n",
        "    - **Intercept** and **slope** describe a **linear** (\"straight line\") relationship between **outcome** $Y_i$ and **predictor variable** $x_i$\n",
        "\n",
        "    \n",
        "- **Error** $\\epsilon_i$ (also sometimes called the **noise**) makes **Simple Linear Regression** a **statistical model** by introducing a **random variable** with a **distribution**\n",
        "\n",
        "    - The $\\sigma^2$ **parameter** is a part of the **noise distribution** and controls how much vertical variability/spread there is in the $Y_i$ data away from the line: $\\sigma^2$ (and $\\beta_0$ as well to be honest) is an \"auxiliary\" **parameter** in the sense that interest is usually in $\\beta_1$ rather than $\\sigma^2$\n",
        "    - **Errors** $\\epsilon_i$ (in conjuction with the **linear form**) define the **assumptions** of the **Simple Linear regression** Model specification\n",
        "    - <u>but these **assumptions** are not the focus of further detailed reviewed here</u>\n",
        "    \n",
        "</details>\n",
        "\n",
        "<details class=\"details-example\"><summary style=\"color:blue\"><u>The implication of the specification on the models assumptions</u></summary>\n",
        "\n",
        "### The implication of the specification on the models assumptions are as follows...\n",
        "    \n",
        "> The first three assumptions associated with the **Simple Linear regression** model are that\n",
        ">\n",
        "> 1. the $\\epsilon_i$ **errors** (sometimes referred to as the **noise**) are **normally distributed**\n",
        "> 2. the $\\epsilon_i$ **errors** are **homoscedastic** (so their distributional variance $\\sigma^2$ does not change as a function of $x_i$)\n",
        "> 3. and the $\\epsilon_i$ **errors** are **statistically independent** (so their values do not depend on each other)\n",
        "> 4. the linear form is [at least reasonably approximately] \"true\" (in the sense that the above two remain [at least reasonably approximately] \"true\") so that then behavior of the $Y_i$ **outcomes** are represented/determined on average by the **linear equation**)<br>\n",
        ">\n",
        ">    and there are just a couple; but, a super deeper reflection on these at this point is \"beyond the scope\" of STA130; nonetheless, they are not too hard to understand and are that<br><br>\n",
        "> 5. the $x_i$ **predictor variable** is assumed to be **measured without error** (so they are taken to have no *randomness* and therefore not have any *distributional* nature: they are just numbers **not random variables**)\n",
        "> 6. and the $\\epsilon_i$ **errors** are **unbiased** relative to the **expected value** of **outcome** $E[Y_i|x_i]=\\beta_0 + \\beta_1x_i$ (which is equivalently stated by saying that the mean of the **error distribution** is $0$, or again equivalently, that the **expected value** of the **errors** $E[\\epsilon_i] = 0$)\n",
        "\n",
        "</details>\n",
        "  \n",
        "<details class=\"details-example\"><summary><u><span style=\"color:blue\">Model fitting for Simple Linear Regression</span></u></summary>\n",
        "\n",
        "### Model fitting for Simple Linear Regression\n",
        "    \n",
        "> Do you remember how to use `statsmodels.ols(...).fit().summary()`?\n",
        ">\n",
        "> Hopefully so, but if not, today's demo with **multiple linear regression** will likely refresh your memory.\n",
        "    \n",
        "The $\\hat y_i = \\hat \\beta_0 + \\hat \\beta_1 x_i$ **fitted model** equation distinctly contrasts with the $Y_i = \\beta_0 + \\beta_1 x_i + \\epsilon_i$ **theoretical model** specification. To emphasize and clarify the difference, we augment our **simple linear regression** model nomenclature (as given in the 'The theoretical \"normal distribution\" regression model' link above) with the contrasting alternative notations and terminology:\n",
        "\n",
        "- **Fitted intercept** $\\hat \\beta_0$ and **slope** $\\hat \\beta_1$ ***coefficients*** are given \"hats\" to distinguish that they **estimate** (based on observed **sample data**), respectively, the **intercept** $\\beta_0$ and **slope** $\\beta_1$ ***parameters***\n",
        "\n",
        "    \n",
        "- **Fitted (predicted) values** $\\hat y_i$ are made lower case and also given \"hats\" to distinguish them from the (upper case) **theoretical random variable** $Y_i$ implied by the **theoretical simple linear regression model**\n",
        "    - Technically, the **error** $\\epsilon_i$ is the **random variable** specified by the **simple linear regression model** specification, and/but this implies the **random variable** nature of $Y_i$\n",
        "\n",
        "    \n",
        "- The **residuals** $\\text{e}_i = \\hat \\epsilon_i = y_i - \\hat y_i = y_i - \\hat \\beta_0 + \\hat \\beta_1 x_i $ are also distinct and definitively contrast with the **errors** (or **noises**) $\\epsilon_i$\n",
        "    - The **residuals** $\\text{e}_i = \\hat \\epsilon_i$ are actually available, while the **error** (or **noises**) $\\epsilon_i$ are just a theoretical concept\n",
        "    - The **residuals** $\\text{e}_i = \\hat \\epsilon_i$ nonetheless are therefore used to diagnostically assess the theoretical modeling assumptions of the  **errors** $\\epsilon_i$, such as the **normality**, **homoskedasticity**, and **linear form** assumptions; and, <u>while this is a not necessarily beyond the scope of STA130 and would certainly be a relevant consideration for the course project, this will not be addressed here at this time</u>\n",
        "\n",
        "</details>\n",
        "      \n",
        "<details class=\"details-example\"><summary><u><span style=\"color:blue\">Hypothesis testing for Simple Linear Regression</span></u></summary>\n",
        "    \n",
        "### Hypothesis testing for Simple Linear Regression\n",
        "    \n",
        "We can use **Simple Linear Regression** to test\n",
        "    \n",
        "$$\\Large\n",
        "\\begin{align}\n",
        "H_0: {}& \\beta_1=0 \\quad \\text{ (there is no linear assocation between $Y_i$ and $x_i$ \"on average\")}\\\\\n",
        "H_A: {}& H_0 \\text{ is false}\n",
        "\\end{align}$$\n",
        "\n",
        "That is, we can assess the evidence of a linear association in the data based on a **null hypothesis** that the **slope** (the \"on average\" change in $Y_i$ per \"single unit\" change in $x_i$) is zero\n",
        "\n",
        "> We are essentially never (or only very rarely in very special circumstances) interested in a **null hypothesis** concerning the **intercept** $\\beta_0$ (as opposed to $\\beta_1$) because the assumption that $\\beta_0$ is zero essentially never (or only very rarely in very special circumstances) has any meaning, whereas the assumption that $\\beta_1$ is zero has the very practically useful interpretation of \"no linear association\" which allows us to evaluate the  evidence of a linear association based on observed data\n",
        "\n",
        "Remember, the **p-value** is \"the probability that a test statistic is as or more extreme than the observed test statistic if the null hypothesis is true\"\n",
        "\n",
        "- We do not prove $H_0$ false, we instead give evidence against the $H_0$\n",
        "     - \"We reject the null hypothesis with a p-value of abc, meaning we have xyz evidence against the null hypothesis\"\n",
        "- We do not prove $H_0$ is true, we instead do not have evidence to reject $H_0$\n",
        "     - \"We fail to reject the null hypothesis with a p-value of abc\"\n",
        "|p-value|Evidence|\n",
        "|-|-|\n",
        "|$$p > 0.1$$|No evidence against the null hypothesis|\n",
        "|$$0.1 \\ge p > 0.05$$|Weak evidence against the null hypothesis|\n",
        "|$$0.05 \\ge p > 0.01$$|Moderate evidence against the null hypothesis|\n",
        "|$$0.01 \\ge p > 0.001$$|Strong evidence against the null hypothesis|\n",
        "|$$0.001 \\ge p$$|Very strong evidence against the null hypothesis|\n",
        "\n",
        "</details>\n",
        "\n",
        "<details class=\"details-example\"><summary><u><span style=\"color:blue\">Indicator variables and (two-sample) group comparision</span></u></summary>\n",
        "\n",
        "### Indicator variables and (two-sample) group comparision\n",
        "    \n",
        "A **Simple Linear Regression** can specify a model of two normally distributed populations (say, \"A\" and \"B\") with different **means** (but a common **variance**)\n",
        "    \n",
        "$$\\Large Y_i = \\beta_0 + \\beta_1 1_{[x_i=\\textrm{\"B\"}]}(x_i) + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "and where the **indicator variable** notation $1_{[x_i=\\textrm{\"B\"}]}(x_i)$ equates to $1$ whenever $x_i$ takes on the value \"B\" (as opposed to population \"A\" in this example), and $0$ otherwise. Thus, the mean of population \"A\" is $\\mu_A = \\beta_0$ while the mean of population \"B\" is $\\mu_A = \\beta_0+\\beta_1$. The $\\beta_1$ coefficient therefore now represents a so-called **contrast** which captures the difference between the two groups (of the \"A\" and \"B\" populations in this example).\n",
        "    \n",
        "A two-group comparison for **independent** (as opposed to **paired**) was also considered from a simulation perspective.\n",
        "\n",
        "- For **simulation-based hypothesis testing** we considered a **permutation test**, which approximated **the sampling distribution for an average difference statistic under the null hypothesis assumption of \"no difference between treatment groups\"**. In **permutation testing** this is done by recalculating the average difference between \"samples\" created by shuffling (or randomly permuting) the label assignments. Shuffling imposes the constraint that the created \"samples\" have the same size as the original two groups. It also fully conforms with the a null hypothesis assumption of \"no difference between treatment groups\" since if this is true the labels don't matter.\n",
        "\n",
        "- A **bootstrapped confidence interval** for the **difference in population means** can also be created by repeatedly resampling each of the samples (separately within each sample) to build up the **bootstrapped sampling distribution of the average difference statistics**.\n",
        "    \n",
        "Students may review the code below in their personal time at their convenience should they wish to review an explicit implementation of the above two-sample analysis specifications.\n",
        "\n",
        "```python\n",
        "# Slight editing of https://chatgpt.com/share/6b1bb97b-80b9-4a50-9323-56a51060d6c8\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "np.random.seed(42)  # Seed for reproducibility\n",
        "# Sample data... more work needed to \"parameterize\" this code to explore different sample sizes, etc.\n",
        "data = {'group': ['A'] * 10 + ['B'] * 10,  # 10 observations for each group... change loc to change group offset\n",
        "        'value': np.random.normal(loc=50, scale=5, size=10).tolist() + np.random.normal(loc=55, scale=5, size=10).tolist()}\n",
        "df = pd.DataFrame(data)\n",
        "\n",
        "# Observed difference in means between group A and B\n",
        "observed_diff = df[df['group'] == 'A']['value'].mean() - df[df['group'] == 'B']['value'].mean()\n",
        "\n",
        "n_permutations = 1000  # Number of permutations\n",
        "permutation_diffs = np.zeros(1000)  # Store the permutation results\n",
        "\n",
        "# Permutation test\n",
        "for i in range(n_permutations):\n",
        "    # Shuffle group labels randomly; Calculate the difference in means with shuffled labels\n",
        "    shuffled = df['group'].sample(frac=1, replace=False).reset_index(drop=True)\n",
        "    shuffled_diff = df[shuffled == 'A']['value'].mean() - df[shuffled == 'B']['value'].mean()\n",
        "    permutation_diffs[i] = shuffled_diff\n",
        "\n",
        "# Calculate p-value: the proportion of permutations with an\n",
        "# absolute difference greater than or equal to the observed difference\n",
        "p_value = np.mean(np.abs(permutation_diffs) >= np.abs(observed_diff))\n",
        "print(f\"Observed difference in means: {observed_diff:.3f}\")\n",
        "print(f\"P-value from permutation test: {p_value:.3f}\")\n",
        "```\n",
        "    \n",
        "```python  \n",
        "import plotly.express as px\n",
        "import plotly.graph_objects as go\n",
        "\n",
        "# Plotting the distribution of permutation differences\n",
        "fig = px.histogram(permutation_diffs, nbins=30, title='Permutation Test: Sampling Distribution of the Difference in Means')\n",
        "fig.update_traces(marker=dict(color=\"lightblue\"), opacity=0.75)\n",
        "\n",
        "# Add a vertical line for the observed difference\n",
        "fig.add_vline(x=observed_diff, line_dash=\"dash\", line_color=\"red\", annotation_text=\"Observed Diff\", annotation_position=\"top right\")\n",
        "\n",
        "# Highlight the area for p-value computation (two-tailed test)\n",
        "fig.add_trace(go.Scatter(\n",
        "    x=permutation_diffs[permutation_diffs >= np.abs(observed_diff)],\n",
        "    y=[0] * len(permutation_diffs[permutation_diffs >= np.abs(observed_diff)]),\n",
        "    mode='markers', marker=dict(color='red'), name=f'Extreme Values (p-value = {p_value:.3f})'))\n",
        "\n",
        "fig.add_trace(go.Scatter(\n",
        "    x=permutation_diffs[permutation_diffs <= -np.abs(observed_diff)],\n",
        "    y=[0] * len(permutation_diffs[permutation_diffs <= -np.abs(observed_diff)]),\n",
        "    mode='markers', marker=dict(color='red'), name='Extreme Values'))\n",
        "\n",
        "# Update layout to make the plot more informative\n",
        "fig.update_layout(xaxis_title='Difference in Means', yaxis_title='Frequency',\n",
        "                  showlegend=True, legend=dict(title=None))\n",
        "fig.show() # USE `fig.show(renderer=\"png\")` FOR ALL GitHub and MarkUs SUBMISSIONS    \n",
        "```\n",
        "\n",
        "```python\n",
        "n_bootstrap = 1000  # Number of bootstrap samples\n",
        "bootstrap_diffs = np.zeros(n_bootstrap)  # Pre-allocate a numpy array to store the bootstrap differences\n",
        "\n",
        "# Bootstrapping to create the sampling distribution of the difference in means\n",
        "for i in range(n_bootstrap):\n",
        "    # Resample with replacement for both groups A and B\n",
        "    boot_A = df[df['group'] == 'A']['value'].sample(frac=1, replace=True).reset_index(drop=True)\n",
        "    boot_B = df[df['group'] == 'B']['value'].sample(frac=1, replace=True).reset_index(drop=True)\n",
        "    \n",
        "    # Calculate the difference in means for the bootstrap sample and store in pre-allocated array\n",
        "    bootstrap_diffs[i] = boot_A.mean() - boot_B.mean()\n",
        "\n",
        "# Calculate the 95% confidence interval (2.5th and 97.5th percentiles)\n",
        "lower_bound = np.percentile(bootstrap_diffs, 2.5)\n",
        "upper_bound = np.percentile(bootstrap_diffs, 97.5)\n",
        "\n",
        "# Print observed difference and confidence interval\n",
        "print(f\"Observed difference in means: {observed_diff:.3f}\")\n",
        "print(f\"95% Confidence Interval: [{lower_bound:.3f}, {upper_bound:.3f}]\")\n",
        "\n",
        "# Plotting the bootstrap distribution\n",
        "fig = px.histogram(bootstrap_diffs, nbins=30, title='Bootstrapped Sampling Distribution of the Difference in Means')\n",
        "fig.update_traces(marker=dict(color=\"lightblue\"), opacity=0.75)\n",
        "\n",
        "# Add lines for the observed difference and confidence interval\n",
        "fig.add_vline(x=observed_diff, line_dash=\"dash\", line_color=\"red\", annotation_text=\"Observed Diff\", annotation_position=\"top right\")\n",
        "fig.add_vline(x=lower_bound, line_dash=\"dot\", line_color=\"green\", annotation_text=\"Lower 95% CI\", annotation_position=\"bottom left\")\n",
        "fig.add_vline(x=upper_bound, line_dash=\"dot\", line_color=\"green\", annotation_text=\"Upper 95% CI\", annotation_position=\"bottom right\")\n",
        "\n",
        "# Update layout to make the plot more informative\n",
        "fig.update_layout(xaxis_title='Difference in Means', yaxis_title='Frequency', showlegend=False)\n",
        "fig.show() # USE `fig.show(renderer=\"png\")` FOR ALL GitHub and MarkUs SUBMISSIONS\n",
        "```    \n",
        "    \n",
        "</details>"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a877d030",
      "metadata": {
        "id": "a877d030"
      },
      "source": [
        "## 💬 🗣️ Communication Activity #1 [35 minutes]\n",
        "\n",
        "Your TA will group you into <u>**SIX**</u> groups, potentially in groups that correspond to **Course Project teams**<br><br>\n",
        "\n",
        "#### *[25 of the 35 minutes]* Consider data of the form given below, with two columns of **continuous** (numeric decimal) data and one column of categorical data (here with at least three different levels), and answer the following questions within your groups\n",
        "\n",
        "|i | study_hours |class_section |exam_score |\n",
        "|:---|:----|:----|:----|\n",
        "|0 |10.9934280 |A |86.530831 |\n",
        "|1 |9.7234711 |A |84.632809|\n",
        "|2 |11.2953770 |B |87.036506 |\n",
        "|3 |13.0460600 |C |97.952866 |\n",
        "|4 |9.5316930 |C |79.749848|\n",
        "\n",
        "#### *[10 of the 35 minutes]* TAs should interrupt the group discussion every few finutes to see if all groups think they have the answer (or want a little more time), and should facilitate group progress by having one of the groups (ideally not the same group over and over) explain the answer to all the other groups (to confirm everyone has the right answer, or help move groups that are stuck along)\n",
        "\n",
        "\n",
        "1. How could you use ONLY TWO **binary indicator variables** in combination to represent the ALL THREE levels (A, B, and C) in the example above?<br><br>\n",
        "\n",
        "    1. Hint 1: If $x_i$ is the `class_section` of **observation** $i$,  what do we know if $1_{[x_i=\\textrm{\"B\"}]}(x_i)$ takes on the value of $1$ (rather than $0$)\n",
        "    2. Hint 2: If $x_i$ is the `class_section` of **observation** $i$,  what do we know if $1_{[x_i=\\textrm{\"C\"}]}(x_i)$ takes on the value of $1$ (rather than $0$)\n",
        "    3. Hint 3: Using THREE **binary indicator variables** would be unnecessarily redundant for encoding the information about which of three groups you were in... can you see why? Suppose both of the above are $0$... what does this situation indicate?<br><br>\n",
        "    \n",
        "2. What are the **means** of the different `class_section` groups in terms of the parameters of the following model specification?<br><br>\n",
        "   $$Y_i = \\beta_0 + 1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_1 + 1_{[x_i=\\textrm{\"C\"}]}(x_i)\\beta_2 + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$<br>\n",
        "\n",
        "3. What is the nature of the data generated under the following model specification if $Y_i$ is the `exam_score` of **observation** $i$, $z_i$ is the value of `study_hours` for **observation** $i$, and $x_i$ is as described above?<br><br>\n",
        "   $$Y_i = \\beta_0 + 1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_1 + 1_{[x_i=\\textrm{\"C\"}]}(x_i)\\beta_2 + \\beta_3 z_i + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$<br>\n",
        "\n",
        "    1. Hint: the model specification below would be described as \"a stright line relationship between the `exam_score` variable and `study_hours` variable with observed data noisily distributed around the line\". So what changes if the indicator variables are included?<br>\n",
        "    \n",
        "    $$Y_i = \\beta_0 + \\beta_3 z_i + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$<br>\n",
        "\n",
        "4. What is the practical interpretation of how `exam_score` changes relative to `class_section` according to the model specification of the previous question if $\\beta_1$ and $\\beta_2$ are not $0$?<br><br>  \n",
        "\n",
        "    1. Hint: what is the meaning if $\\beta_1$ and $\\beta_2$ are actually $0$?<br><br>  \n",
        "    \n",
        "5. What is the practical interpretation of the behavior of the relationship between `exam_score` and `study_hours` within different `class_section` groups according to the model specification of the previous question?<br><br>\n",
        "\n",
        "    1. Hint: Does the nature of the relationship really change? Does the model specification prescribe that \"the change in the outcome variable on average per unit change in the predictor\" should differ across different `class_section` groups?<br><br>  \n",
        "\n",
        "6. Is there a different kind of behavior that could be seen for the relationship between `exam_score` and `study_hours` between different `class_section` groups that might be different than what's prescribed by the model specification of the previous question?<br><br>  \n",
        "\n",
        "    1. Hint 1: what is the meaning of the following model specification?<br>\n",
        "    \n",
        "    $$Y_i = \\beta_0 + \\beta_3 z_i + 1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_1 + \\beta_4 z_i \\times 1_{[x_i=\\textrm{\"B\"}]}(x_i) + 1_{[x_i=\\textrm{\"C\"}]}(x_i)\\beta_2 + \\beta_5 z_i \\times 1_{[x_i=\\textrm{\"C\"}]}(x_i) + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$<br>\n",
        "\n",
        "    1. Hint 2: this could also be re-expressed as...<br>\n",
        "    \n",
        "    $$Y_i = \\left(\\beta_0 +  1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_1 + 1_{[x_i=\\textrm{\"C\"}]}(x_i)\\beta_2\\right) + \\left(\\beta_3 + \\beta_4 1_{[x_i=\\textrm{\"B\"}]}(x_i) + \\beta_5 \\right) \\times z_i + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import statsmodels.api as sm\n",
        "import statsmodels.formula.api as smf\n",
        "\n",
        "# Given data\n",
        "data = {\n",
        "    \"study_hours\": [10.9934280, 9.7234711, 11.2953770, 13.0460600, 9.5316930],\n",
        "    \"class_section\": [\"A\", \"A\", \"B\", \"C\", \"C\"],\n",
        "    \"exam_score\": [86.530831, 84.632809, 87.036506, 97.952866, 79.749848]\n",
        "}\n",
        "\n",
        "# Create a DataFrame\n",
        "df = pd.DataFrame(data)\n",
        "\n",
        "# Step 1: Create indicator variables for class_section\n",
        "df['is_B'] = (df['class_section'] == 'B').astype(int)\n",
        "df['is_C'] = (df['class_section'] == 'C').astype(int)\n",
        "\n",
        "# Step 2: Fit the basic model without interaction terms\n",
        "model_basic = smf.ols('exam_score ~ is_B + is_C + study_hours', data=df).fit()\n",
        "\n",
        "# Step 3: Fit the model with interaction terms\n",
        "model_interaction = smf.ols('exam_score ~ is_B * study_hours + is_C * study_hours', data=df).fit()\n",
        "\n",
        "# Display the summaries of both models\n",
        "basic_summary = model_basic.summary()\n",
        "interaction_summary = model_interaction.summary()\n",
        "\n",
        "basic_summary, interaction_summary\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "G2Oouu60zZNC",
        "outputId": "4b85ce45-25d2-44ce-945d-07ed811d8e11"
      },
      "id": "G2Oouu60zZNC",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.\n",
            "  warn(\"omni_normtest is not valid with less than 8 observations; %i \"\n",
            "/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.\n",
            "  warn(\"omni_normtest is not valid with less than 8 observations; %i \"\n",
            "/usr/local/lib/python3.10/dist-packages/statsmodels/regression/linear_model.py:1795: RuntimeWarning: divide by zero encountered in divide\n",
            "  return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid)\n",
            "/usr/local/lib/python3.10/dist-packages/statsmodels/regression/linear_model.py:1795: RuntimeWarning: invalid value encountered in scalar multiply\n",
            "  return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid)\n",
            "/usr/local/lib/python3.10/dist-packages/statsmodels/regression/linear_model.py:1717: RuntimeWarning: divide by zero encountered in scalar divide\n",
            "  return np.dot(wresid, wresid) / self.df_resid\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(<class 'statsmodels.iolib.summary.Summary'>\n",
              " \"\"\"\n",
              "                             OLS Regression Results                            \n",
              " ==============================================================================\n",
              " Dep. Variable:             exam_score   R-squared:                       0.946\n",
              " Model:                            OLS   Adj. R-squared:                  0.783\n",
              " Method:                 Least Squares   F-statistic:                     5.799\n",
              " Date:                Fri, 08 Nov 2024   Prob (F-statistic):              0.294\n",
              " Time:                        17:23:57   Log-Likelihood:                -8.7477\n",
              " No. Observations:                   5   AIC:                             25.50\n",
              " Df Residuals:                       1   BIC:                             23.93\n",
              " Df Model:                           3                                         \n",
              " Covariance Type:            nonrobust                                         \n",
              " ===============================================================================\n",
              "                   coef    std err          t      P>|t|      [0.025      0.975]\n",
              " -------------------------------------------------------------------------------\n",
              " Intercept      36.3380     12.397      2.931      0.209    -121.186     193.862\n",
              " is_B           -2.9994      3.968     -0.756      0.588     -53.420      47.421\n",
              " is_C           -1.1537      3.299     -0.350      0.786     -43.078      40.770\n",
              " study_hours     4.7540      1.178      4.036      0.155     -10.212      19.720\n",
              " ==============================================================================\n",
              " Omnibus:                          nan   Durbin-Watson:                   2.500\n",
              " Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.213\n",
              " Skew:                           0.000   Prob(JB):                        0.899\n",
              " Kurtosis:                       1.989   Cond. No.                         99.2\n",
              " ==============================================================================\n",
              " \n",
              " Notes:\n",
              " [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
              " \"\"\",\n",
              " <class 'statsmodels.iolib.summary.Summary'>\n",
              " \"\"\"\n",
              "                             OLS Regression Results                            \n",
              " ==============================================================================\n",
              " Dep. Variable:             exam_score   R-squared:                       1.000\n",
              " Model:                            OLS   Adj. R-squared:                    nan\n",
              " Method:                 Least Squares   F-statistic:                       nan\n",
              " Date:                Fri, 08 Nov 2024   Prob (F-statistic):                nan\n",
              " Time:                        17:23:57   Log-Likelihood:                 140.70\n",
              " No. Observations:                   5   AIC:                            -271.4\n",
              " Df Residuals:                       0   BIC:                            -273.4\n",
              " Df Model:                           4                                         \n",
              " Covariance Type:            nonrobust                                         \n",
              " ====================================================================================\n",
              "                        coef    std err          t      P>|t|      [0.025      0.975]\n",
              " ------------------------------------------------------------------------------------\n",
              " Intercept           70.1005        inf          0        nan         nan         nan\n",
              " is_B                 0.0004        inf          0        nan         nan         nan\n",
              " study_hours          1.4946        inf          0        nan         nan         nan\n",
              " is_B:study_hours     0.0048        inf          0        nan         nan         nan\n",
              " is_C               -39.7211        inf         -0        nan         nan         nan\n",
              " is_C:study_hours     3.6850        inf          0        nan         nan         nan\n",
              " ==============================================================================\n",
              " Omnibus:                          nan   Durbin-Watson:                   0.855\n",
              " Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.709\n",
              " Skew:                           0.911   Prob(JB):                        0.701\n",
              " Kurtosis:                       2.707   Cond. No.                         465.\n",
              " ==============================================================================\n",
              " \n",
              " Notes:\n",
              " [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
              " [2] The input rank is higher than the number of observations.\n",
              " \"\"\")"
            ]
          },
          "metadata": {},
          "execution_count": 1
        }
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6d2c0b98",
      "metadata": {
        "id": "6d2c0b98"
      },
      "source": [
        "##  🚧 🏗️ Demo (combining different kinds of predictor variables) [40 minutes]\n",
        "\n",
        "### Demo statsmodels multiple linear regression \"formula\" specifications [16 of the 40 minutes]\n",
        "\n",
        "- Get a dataset and briefly show the `statsmodels` \"formula\" implementions (https://www.statsmodels.org/dev/example_formulas.html) of the specifications below as you progress through introducing and discussing them; so,\n",
        "    - introduce, explain, discuss, and then `statsmodels` \"formula\" demo them quickly\n",
        "    - (and use a ChatBot to figure things out whenever needed)\n",
        "    \n",
        "### Explain notation, terminology, and meaning [24 of the 40 minutes]\n",
        "\n",
        "- **_There won't be time to figure out the visualizations on the fly so TAs won't show these using `plotly` or `python` (unless they've pre-populated such demos); BUT, the TAs might COULD MAYBE just draw the illustrative pictures on the board?_**\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# Load the dataset\n",
        "file_path = \"/content/sample_data/california_housing_test.csv\"\n",
        "df = pd.read_csv(file_path)\n",
        "\n",
        "# Take the first 200 rows\n",
        "df_subset = df.head(200)\n",
        "\n",
        "# Create a new categorical column with random values 'A', 'B', or 'C'\n",
        "df_subset['category'] = np.random.choice(['A', 'B', 'C'], size=len(df_subset))\n",
        "\n",
        "# Display the first few rows of the new DataFrame\n",
        "print(df_subset.head())\n",
        "\n",
        "# Optionally, save the new DataFrame to a CSV file\n",
        "df_subset.to_csv(\"subset_with_categories.csv\", index=False)\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "h5BrX15e7C2G",
        "outputId": "35022972-04ee-4f34-87b4-d84eeb6301e9"
      },
      "id": "h5BrX15e7C2G",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "   longitude  latitude  housing_median_age  total_rooms  total_bedrooms  \\\n",
            "0    -122.05     37.37                27.0       3885.0           661.0   \n",
            "1    -118.30     34.26                43.0       1510.0           310.0   \n",
            "2    -117.81     33.78                27.0       3589.0           507.0   \n",
            "3    -118.36     33.82                28.0         67.0            15.0   \n",
            "4    -119.67     36.33                19.0       1241.0           244.0   \n",
            "\n",
            "   population  households  median_income  median_house_value category  \n",
            "0      1537.0       606.0         6.6085            344700.0        C  \n",
            "1       809.0       277.0         3.5990            176500.0        A  \n",
            "2      1484.0       495.0         5.7934            270500.0        A  \n",
            "3        49.0        11.0         6.1359            330000.0        C  \n",
            "4       850.0       237.0         2.9375             81700.0        C  \n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-7-81ac2a27e541>:12: SettingWithCopyWarning: \n",
            "A value is trying to be set on a copy of a slice from a DataFrame.\n",
            "Try using .loc[row_indexer,col_indexer] = value instead\n",
            "\n",
            "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
            "  df_subset['category'] = np.random.choice(['A', 'B', 'C'], size=len(df_subset))\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "file_path = \"/content/subset_with_categories.csv\"\n",
        "df = pd.read_csv(file_path)\n",
        "df.head(200)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 444
        },
        "id": "OLFeWuUS9_PF",
        "outputId": "c8ce0984-50dc-46bd-f036-812ebf74bd70"
      },
      "id": "OLFeWuUS9_PF",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "     longitude  latitude  housing_median_age  total_rooms  total_bedrooms  \\\n",
              "0      -122.05     37.37                27.0       3885.0           661.0   \n",
              "1      -118.30     34.26                43.0       1510.0           310.0   \n",
              "2      -117.81     33.78                27.0       3589.0           507.0   \n",
              "3      -118.36     33.82                28.0         67.0            15.0   \n",
              "4      -119.67     36.33                19.0       1241.0           244.0   \n",
              "..         ...       ...                 ...          ...             ...   \n",
              "195    -116.92     32.86                11.0       2204.0           518.0   \n",
              "196    -121.04     38.95                22.0       1931.0           445.0   \n",
              "197    -122.12     37.45                38.0       1276.0           314.0   \n",
              "198    -119.48     36.54                28.0       2112.0           363.0   \n",
              "199    -121.02     37.68                25.0       3262.0           588.0   \n",
              "\n",
              "     population  households  median_income  median_house_value category  \n",
              "0        1537.0       606.0         6.6085            344700.0        C  \n",
              "1         809.0       277.0         3.5990            176500.0        A  \n",
              "2        1484.0       495.0         5.7934            270500.0        A  \n",
              "3          49.0        11.0         6.1359            330000.0        C  \n",
              "4         850.0       237.0         2.9375             81700.0        C  \n",
              "..          ...         ...            ...                 ...      ...  \n",
              "195      1472.0       497.0         2.3693            127000.0        A  \n",
              "196      1009.0       407.0         2.7500            153200.0        A  \n",
              "197       955.0       287.0         2.0096            155700.0        C  \n",
              "198      1011.0       335.0         4.2222            108900.0        B  \n",
              "199      1834.0       578.0         3.9960            114500.0        A  \n",
              "\n",
              "[200 rows x 10 columns]"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-52dc2d15-dfd5-4c14-8c02-cefd78752426\" class=\"colab-df-container\">\n",
              "    <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>longitude</th>\n",
              "      <th>latitude</th>\n",
              "      <th>housing_median_age</th>\n",
              "      <th>total_rooms</th>\n",
              "      <th>total_bedrooms</th>\n",
              "      <th>population</th>\n",
              "      <th>households</th>\n",
              "      <th>median_income</th>\n",
              "      <th>median_house_value</th>\n",
              "      <th>category</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>-122.05</td>\n",
              "      <td>37.37</td>\n",
              "      <td>27.0</td>\n",
              "      <td>3885.0</td>\n",
              "      <td>661.0</td>\n",
              "      <td>1537.0</td>\n",
              "      <td>606.0</td>\n",
              "      <td>6.6085</td>\n",
              "      <td>344700.0</td>\n",
              "      <td>C</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>-118.30</td>\n",
              "      <td>34.26</td>\n",
              "      <td>43.0</td>\n",
              "      <td>1510.0</td>\n",
              "      <td>310.0</td>\n",
              "      <td>809.0</td>\n",
              "      <td>277.0</td>\n",
              "      <td>3.5990</td>\n",
              "      <td>176500.0</td>\n",
              "      <td>A</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>-117.81</td>\n",
              "      <td>33.78</td>\n",
              "      <td>27.0</td>\n",
              "      <td>3589.0</td>\n",
              "      <td>507.0</td>\n",
              "      <td>1484.0</td>\n",
              "      <td>495.0</td>\n",
              "      <td>5.7934</td>\n",
              "      <td>270500.0</td>\n",
              "      <td>A</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>-118.36</td>\n",
              "      <td>33.82</td>\n",
              "      <td>28.0</td>\n",
              "      <td>67.0</td>\n",
              "      <td>15.0</td>\n",
              "      <td>49.0</td>\n",
              "      <td>11.0</td>\n",
              "      <td>6.1359</td>\n",
              "      <td>330000.0</td>\n",
              "      <td>C</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>-119.67</td>\n",
              "      <td>36.33</td>\n",
              "      <td>19.0</td>\n",
              "      <td>1241.0</td>\n",
              "      <td>244.0</td>\n",
              "      <td>850.0</td>\n",
              "      <td>237.0</td>\n",
              "      <td>2.9375</td>\n",
              "      <td>81700.0</td>\n",
              "      <td>C</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>195</th>\n",
              "      <td>-116.92</td>\n",
              "      <td>32.86</td>\n",
              "      <td>11.0</td>\n",
              "      <td>2204.0</td>\n",
              "      <td>518.0</td>\n",
              "      <td>1472.0</td>\n",
              "      <td>497.0</td>\n",
              "      <td>2.3693</td>\n",
              "      <td>127000.0</td>\n",
              "      <td>A</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>196</th>\n",
              "      <td>-121.04</td>\n",
              "      <td>38.95</td>\n",
              "      <td>22.0</td>\n",
              "      <td>1931.0</td>\n",
              "      <td>445.0</td>\n",
              "      <td>1009.0</td>\n",
              "      <td>407.0</td>\n",
              "      <td>2.7500</td>\n",
              "      <td>153200.0</td>\n",
              "      <td>A</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>197</th>\n",
              "      <td>-122.12</td>\n",
              "      <td>37.45</td>\n",
              "      <td>38.0</td>\n",
              "      <td>1276.0</td>\n",
              "      <td>314.0</td>\n",
              "      <td>955.0</td>\n",
              "      <td>287.0</td>\n",
              "      <td>2.0096</td>\n",
              "      <td>155700.0</td>\n",
              "      <td>C</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>198</th>\n",
              "      <td>-119.48</td>\n",
              "      <td>36.54</td>\n",
              "      <td>28.0</td>\n",
              "      <td>2112.0</td>\n",
              "      <td>363.0</td>\n",
              "      <td>1011.0</td>\n",
              "      <td>335.0</td>\n",
              "      <td>4.2222</td>\n",
              "      <td>108900.0</td>\n",
              "      <td>B</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>199</th>\n",
              "      <td>-121.02</td>\n",
              "      <td>37.68</td>\n",
              "      <td>25.0</td>\n",
              "      <td>3262.0</td>\n",
              "      <td>588.0</td>\n",
              "      <td>1834.0</td>\n",
              "      <td>578.0</td>\n",
              "      <td>3.9960</td>\n",
              "      <td>114500.0</td>\n",
              "      <td>A</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>200 rows × 10 columns</p>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-52dc2d15-dfd5-4c14-8c02-cefd78752426')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-52dc2d15-dfd5-4c14-8c02-cefd78752426 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-52dc2d15-dfd5-4c14-8c02-cefd78752426');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-164caf8d-ed64-4546-8150-247a420fa3b0\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-164caf8d-ed64-4546-8150-247a420fa3b0')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-164caf8d-ed64-4546-8150-247a420fa3b0 button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "dataframe",
              "variable_name": "df",
              "summary": "{\n  \"name\": \"df\",\n  \"rows\": 200,\n  \"fields\": [\n    {\n      \"column\": \"longitude\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.9277857321667606,\n        \"min\": -124.16,\n        \"max\": -116.92,\n        \"num_unique_values\": 157,\n        \"samples\": [\n          -121.64,\n          -118.29,\n          -117.7\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"latitude\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2.106230751220516,\n        \"min\": 32.72,\n        \"max\": 41.01,\n        \"num_unique_values\": 137,\n        \"samples\": [\n          36.95,\n          40.46,\n          37.25\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"housing_median_age\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 11.730344406687193,\n        \"min\": 3.0,\n        \"max\": 52.0,\n        \"num_unique_values\": 47,\n        \"samples\": [\n          4.0,\n          50.0,\n          48.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"total_rooms\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1960.683576308792,\n        \"min\": 67.0,\n        \"max\": 16181.0,\n        \"num_unique_values\": 199,\n        \"samples\": [\n          1583.0,\n          161.0,\n          2489.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"total_bedrooms\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 371.1178956035211,\n        \"min\": 15.0,\n        \"max\": 2971.0,\n        \"num_unique_values\": 181,\n        \"samples\": [\n          1307.0,\n          226.0,\n          307.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"population\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1088.8693692502588,\n        \"min\": 49.0,\n        \"max\": 8152.0,\n        \"num_unique_values\": 194,\n        \"samples\": [\n          1645.0,\n          189.0,\n          2380.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"households\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 338.16245148295934,\n        \"min\": 11.0,\n        \"max\": 2651.0,\n        \"num_unique_values\": 172,\n        \"samples\": [\n          351.0,\n          663.0,\n          362.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"median_income\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 2.2132031635481852,\n        \"min\": 0.4999,\n        \"max\": 15.0001,\n        \"num_unique_values\": 197,\n        \"samples\": [\n          3.7727,\n          2.1422,\n          1.875\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"median_house_value\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 105905.38704790993,\n        \"min\": 47500.0,\n        \"max\": 500001.0,\n        \"num_unique_values\": 177,\n        \"samples\": [\n          126300.0,\n          212500.0,\n          112500.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"category\",\n      \"properties\": {\n        \"dtype\": \"category\",\n        \"num_unique_values\": 3,\n        \"samples\": [\n          \"C\",\n          \"A\",\n          \"B\"\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}"
            }
          },
          "metadata": {},
          "execution_count": 9
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c25a0cef",
      "metadata": {
        "id": "c25a0cef"
      },
      "outputs": [],
      "source": [
        "# Hopefully you're TA picks an interesting data set to demo!\n",
        "# Students *could* potentially help by suggesting an intersting one to try out, too..."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "93921cd1",
      "metadata": {
        "id": "93921cd1"
      },
      "source": [
        "1. Two (or more) **indicator predictor variables** corresponding to THREE (or more) groups based on so called \"contrasts\" (or \"offsets\" or differnces) from an (arbitrarily chosen) \"baseline\" group\n",
        "\n",
        "   $$\\large Y_i = \\beta_{\\textrm{A}} + 1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_{\\textrm{B-offset}} + 1_{[x_i=\\textrm{\"C\"}]}(x_i)\\beta_{\\textrm{C-offset}} + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import statsmodels.api as sm\n",
        "import statsmodels.formula.api as smf\n",
        "\n",
        "import numpy as np\n",
        "\n",
        "import pandas\n"
      ],
      "metadata": {
        "id": "NAKP_0mFfHRm"
      },
      "id": "NAKP_0mFfHRm",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ee4df5c4",
      "metadata": {
        "id": "ee4df5c4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e2d79030-a0a6-4a70-aad0-64e5df7a5d09"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "                            OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:     median_house_value   R-squared:                       0.003\n",
            "Model:                            OLS   Adj. R-squared:                 -0.002\n",
            "Method:                 Least Squares   F-statistic:                    0.6468\n",
            "Date:                Fri, 08 Nov 2024   Prob (F-statistic):              0.422\n",
            "Time:                        20:35:03   Log-Likelihood:                -2597.0\n",
            "No. Observations:                 200   AIC:                             5198.\n",
            "Df Residuals:                     198   BIC:                             5205.\n",
            "Df Model:                           1                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "==================================================================================\n",
            "                     coef    std err          t      P>|t|      [0.025      0.975]\n",
            "----------------------------------------------------------------------------------\n",
            "Intercept       1.888e+05   1.29e+04     14.641      0.000    1.63e+05    2.14e+05\n",
            "total_bedrooms    16.2835     20.247      0.804      0.422     -23.644      56.211\n",
            "==============================================================================\n",
            "Omnibus:                       37.512   Durbin-Watson:                   1.809\n",
            "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               52.302\n",
            "Skew:                           1.146   Prob(JB):                     4.39e-12\n",
            "Kurtosis:                       4.009   Cond. No.                     1.10e+03\n",
            "==============================================================================\n",
            "\n",
            "Notes:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "[2] The condition number is large, 1.1e+03. This might indicate that there are\n",
            "strong multicollinearity or other numerical problems.\n"
          ]
        }
      ],
      "source": [
        "# Anyway your TA's gonna be live coding for ya here...\n",
        "# continuous variable only: x: total bedrooms  y : median house value\n",
        "mod = smf.ols(formula='median_house_value ~ total_bedrooms', data=df)\n",
        "\n",
        "res = mod.fit()\n",
        "\n",
        "print(res.summary())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "613c9c0f",
      "metadata": {
        "id": "613c9c0f"
      },
      "source": [
        "2. Both **indicator** and **continuous predictor variables** together in a \"parallel lines\" specification where the number of such lines increases by one for each **indicator predictor variables**\n",
        "\n",
        "  $$\\large Y_i = \\beta_{\\textrm{A}} + 1_{[x_i=\\textrm{\"B\"}]}(x_i)\\beta_{\\textrm{B-offset}} + \\beta_z z_{i} + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "    - Within each of the groups induced by the **indicator variables(s)**, the coefficient $\\beta_z$ is (of course) still interpreted as a slope, as \"the average change in the outcome variable for a one-unit increase in the $z_{k}$ predictor variable\" but there is now also an additional difference \"on average\" (of $\\beta_{\\textrm{B-offset}}$ between groups \"A\" and \"B\" in this example) which is what produces the \"parallel lines\" nature of this specification"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a9296c2f",
      "metadata": {
        "id": "a9296c2f",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "be85ff4a-9a4a-4417-b8e0-995b9d35481f"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "                            OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:     median_house_value   R-squared:                       0.007\n",
            "Model:                            OLS   Adj. R-squared:                 -0.009\n",
            "Method:                 Least Squares   F-statistic:                    0.4409\n",
            "Date:                Fri, 08 Nov 2024   Prob (F-statistic):              0.724\n",
            "Time:                        20:37:09   Log-Likelihood:                -2596.7\n",
            "No. Observations:                 200   AIC:                             5201.\n",
            "Df Residuals:                     196   BIC:                             5215.\n",
            "Df Model:                           3                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "====================================================================================\n",
            "                       coef    std err          t      P>|t|      [0.025      0.975]\n",
            "------------------------------------------------------------------------------------\n",
            "Intercept         1.832e+05   1.71e+04     10.715      0.000    1.49e+05    2.17e+05\n",
            "C(category)[T.B]  1318.3243   1.93e+04      0.068      0.946   -3.68e+04    3.95e+04\n",
            "C(category)[T.C]  1.331e+04   1.76e+04      0.754      0.452   -2.15e+04    4.81e+04\n",
            "total_bedrooms      16.7046     20.400      0.819      0.414     -23.527      56.936\n",
            "==============================================================================\n",
            "Omnibus:                       38.752   Durbin-Watson:                   1.816\n",
            "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               54.950\n",
            "Skew:                           1.163   Prob(JB):                     1.17e-12\n",
            "Kurtosis:                       4.087   Cond. No.                     2.20e+03\n",
            "==============================================================================\n",
            "\n",
            "Notes:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "[2] The condition number is large, 2.2e+03. This might indicate that there are\n",
            "strong multicollinearity or other numerical problems.\n"
          ]
        }
      ],
      "source": [
        "# here\n",
        "# continuous variable + categorical variable\n",
        "mod = smf.ols(formula='median_house_value ~ total_bedrooms + C(category)', data=df)\n",
        "\n",
        "res = mod.fit()\n",
        "\n",
        "print(res.summary())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c330770f",
      "metadata": {
        "id": "c330770f"
      },
      "source": [
        "3. Two (or more) **continuous predictor variables**, say $z_{1}$ and $z_{2}$ (with $z_{1i}$ and $z_{2i}$ representing observation $i$), specifying simultaneous linear relationships with the outcome variable occurring in concert\n",
        "\n",
        "  $$\\large Y_i = \\beta_0 + \\beta_1 z_{1i} + \\beta_2 z_{2i} + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "    - Coefficients $\\beta_k$ are still interpreted as \"the average change in the outcome variable for a one-unit increase in the $z_{k}$ predictor variable\" but now we should be careful to additionally add the stipulation \"with all other predictor variables held constant\"\n",
        "\n",
        "|Which actually in fact just happens to be a (hyper)plane model|\n",
        "|-|\n",
        "|![](https://statsandr.com/blog/multiple-linear-regression-made-simple/images/multiple-linear-regression-plane.png)|"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8735c876",
      "metadata": {
        "id": "8735c876",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a61eea4e-72f8-4ab1-e6f1-8d1a23525ebc"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "                            OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:     median_house_value   R-squared:                       0.021\n",
            "Model:                            OLS   Adj. R-squared:                  0.011\n",
            "Method:                 Least Squares   F-statistic:                     2.073\n",
            "Date:                Fri, 08 Nov 2024   Prob (F-statistic):              0.129\n",
            "Time:                        20:40:03   Log-Likelihood:                -2595.3\n",
            "No. Observations:                 200   AIC:                             5197.\n",
            "Df Residuals:                     197   BIC:                             5206.\n",
            "Df Model:                           2                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "==================================================================================\n",
            "                     coef    std err          t      P>|t|      [0.025      0.975]\n",
            "----------------------------------------------------------------------------------\n",
            "Intercept       1.932e+05    1.3e+04     14.830      0.000    1.67e+05    2.19e+05\n",
            "total_bedrooms    84.2435     41.564      2.027      0.044       2.276     166.211\n",
            "population       -26.4712     14.166     -1.869      0.063     -54.408       1.466\n",
            "==============================================================================\n",
            "Omnibus:                       40.440   Durbin-Watson:                   1.850\n",
            "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               59.244\n",
            "Skew:                           1.172   Prob(JB):                     1.37e-13\n",
            "Kurtosis:                       4.271   Cond. No.                     3.40e+03\n",
            "==============================================================================\n",
            "\n",
            "Notes:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "[2] The condition number is large, 3.4e+03. This might indicate that there are\n",
            "strong multicollinearity or other numerical problems.\n"
          ]
        }
      ],
      "source": [
        "# here\n",
        "# 2 continuous variable with no interaction\n",
        "mod = smf.ols(formula='median_house_value ~ total_bedrooms + population\t', data=df)\n",
        "\n",
        "res = mod.fit()\n",
        "\n",
        "print(res.summary())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a83a5ca9",
      "metadata": {
        "id": "a83a5ca9"
      },
      "source": [
        "4. **Interactions** (not to be confused with **indicator predictor variables**) between **continuous** and **indicator predictor variables** allowing the strength of the \"linear relationship\" between the **outcome variable** and the **continuous predictor variables** to be different within different groups, creating a \"non-parallel lines\" behavior\n",
        "\n",
        "    $$\\large Y_i = \\beta_0 + \\beta_z z_i + \\beta_{\\textrm{B-offset}} 1_{[x_i=\\textrm{\"B\"}]}(x_i) + \\beta_{\\textrm{z-change-in-B}} z_i \\times 1_{[x_i=\\textrm{\"B\"}]}(x_i) + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$<br>\n",
        "    $$\\Large\\textrm{OR}$$\n",
        "    $$\\large Y_i = \\beta_0 + \\beta_{\\textrm{B-offset}} 1_{[x_i=\\textrm{\"B\"}]}(x_i) + \\left(\\beta_z +\\beta_{\\textrm{z-change-in-B}} 1_{[x_i=\\textrm{\"B\"}]}(x_i) \\right) \\times z_i + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "    - So the coefficient $\\beta_z$ is (still) interpreted as \"the average change in the outcome variable for a one-unit increase in the $z$ predictor variable within the baseline group\" but now this \"average change\" is different in other groups, so (for this example), \"the average change in the outcome variable for a one-unit increase in the $z$ predictor variable WITHIN GROUP 'B' is $\\beta_z + \\beta_{\\textrm{z-change-in-B}}$\"<br><br>\n",
        "    - And there remains here as well an additional difference \"on average\" (of $\\beta_{\\textrm{B-offset}}$ between groups \"A\" and \"B\" in this example) allowing the \"non-parallel lines\" to have different vertical shifts vertically as well as have different slopes\n",
        "    \n",
        "> When using [`plotly.express.scatter`](https://plotly.com/python/linear-fits/#fitting-multiple-lines-and-retrieving-the-model-parameters) with `x`, `y`, AND `color` AND `trendline=\"ols\"` you automatically get the \"interaction\" visualization."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "af0da2f2",
      "metadata": {
        "id": "af0da2f2",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "b3f2595f-6a59-4338-a7b6-95b57d41ba22"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "                            OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:     median_house_value   R-squared:                       0.028\n",
            "Model:                            OLS   Adj. R-squared:                  0.003\n",
            "Method:                 Least Squares   F-statistic:                     1.100\n",
            "Date:                Fri, 08 Nov 2024   Prob (F-statistic):              0.362\n",
            "Time:                        20:46:20   Log-Likelihood:                -2594.6\n",
            "No. Observations:                 200   AIC:                             5201.\n",
            "Df Residuals:                     194   BIC:                             5221.\n",
            "Df Model:                           5                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "===================================================================================================\n",
            "                                      coef    std err          t      P>|t|      [0.025      0.975]\n",
            "---------------------------------------------------------------------------------------------------\n",
            "Intercept                        1.998e+05   1.97e+04     10.119      0.000    1.61e+05    2.39e+05\n",
            "C(category)[T.B]                -5.889e+04   3.56e+04     -1.655      0.100   -1.29e+05    1.13e+04\n",
            "C(category)[T.C]                -9948.3946   2.99e+04     -0.333      0.740   -6.89e+04    4.91e+04\n",
            "total_bedrooms                    -12.9347     27.066     -0.478      0.633     -66.315      40.446\n",
            "total_bedrooms:C(category)[T.B]   121.3257     60.762      1.997      0.047       1.487     241.165\n",
            "total_bedrooms:C(category)[T.C]    42.6894     45.925      0.930      0.354     -47.887     133.266\n",
            "==============================================================================\n",
            "Omnibus:                       36.958   Durbin-Watson:                   1.881\n",
            "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               51.528\n",
            "Skew:                           1.121   Prob(JB):                     6.47e-12\n",
            "Kurtosis:                       4.076   Cond. No.                     3.85e+03\n",
            "==============================================================================\n",
            "\n",
            "Notes:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "[2] The condition number is large, 3.85e+03. This might indicate that there are\n",
            "strong multicollinearity or other numerical problems.\n"
          ]
        }
      ],
      "source": [
        "# here\n",
        "#  continuous variable with interaction to catogorical variable\n",
        "mod = smf.ols(formula='median_house_value ~ total_bedrooms * C(category)\t', data=df)\n",
        "\n",
        "res = mod.fit()\n",
        "\n",
        "print(res.summary())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "41550929",
      "metadata": {
        "id": "41550929"
      },
      "source": [
        "\n",
        "5. **Interactions** between **continuous predictor variables** creating a \"synergistic\" behavior between the **predictor variables** by which \"the average change in the outcome variable for a one-unit increase in one $z_k$ predictor variable depends on the value of another $z_{k'}$ predictor variable (assuming this other $z_{k'}$ predictor variable is held constant)\"<br><br>\n",
        "  $$\\large Y_i = \\beta_0 + \\beta_1 z_{1i} + \\beta_2 z_{2i} + \\beta_{12} z_{1i} \\times z_{2i} + \\epsilon_i \\quad \\text{ where } \\quad \\epsilon_i \\sim \\mathcal N\\left(0, \\sigma^2\\right)$$\n",
        "\n",
        "    - For this example, \"the average change in the outcome variable for a one-unit increase in the $z_{1}$  predictor variable (assuming $z_{2}$ predictor variable is held constant) is $\\beta_1 + \\beta_{12} z_{2}$\" which can be seen to \"synergistically\" depend on the (fixed) value of $z_{2}$\n",
        "    - And conversely, \"the average change in the outcome variable for a one-unit increase in the $z_{2}$  predictor variable (assuming $z_{1}$ predictor variable is held constant) is $\\beta_1 + \\beta_{12} z_{1}$\" which can be seen to \"synergistically\" depend on the (fixed) value of $z_{1}$<br><br>\n",
        "    - But do you see why?\n",
        "\n",
        "> A quite interesting **continuous predictor variables** \"interaction\" is that of *a variable with itsef*, such as $z_{1}^2$ and $z_{2}^2$.  \n",
        ">\n",
        "> ![](https://miro.medium.com/v2/resize:fit:1200/1*_UaCxPswsCxkj9JzYXCiWg.png)\n",
        "> ![](https://www.researchgate.net/publication/316420362/figure/fig4/AS:486432612589585@1492985980566/First-eight-terms-of-a-2-D-Chebyshev-polynomial.png)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ecad50d8",
      "metadata": {
        "id": "ecad50d8",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "9407ce12-aad9-46f5-c1f2-f9b2f9a4049b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "                            OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:     median_house_value   R-squared:                       0.032\n",
            "Model:                            OLS   Adj. R-squared:                  0.017\n",
            "Method:                 Least Squares   F-statistic:                     2.139\n",
            "Date:                Fri, 08 Nov 2024   Prob (F-statistic):             0.0966\n",
            "Time:                        20:43:10   Log-Likelihood:                -2594.1\n",
            "No. Observations:                 200   AIC:                             5196.\n",
            "Df Residuals:                     196   BIC:                             5209.\n",
            "Df Model:                           3                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "=============================================================================================\n",
            "                                coef    std err          t      P>|t|      [0.025      0.975]\n",
            "---------------------------------------------------------------------------------------------\n",
            "Intercept                  1.721e+05   1.91e+04      8.987      0.000    1.34e+05     2.1e+05\n",
            "total_bedrooms              127.5640     50.526      2.525      0.012      27.919     227.209\n",
            "population                  -18.1301     15.180     -1.194      0.234     -48.067      11.806\n",
            "total_bedrooms:population    -0.0123      0.008     -1.498      0.136      -0.028       0.004\n",
            "==============================================================================\n",
            "Omnibus:                       40.568   Durbin-Watson:                   1.877\n",
            "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               59.607\n",
            "Skew:                           1.172   Prob(JB):                     1.14e-13\n",
            "Kurtosis:                       4.287   Cond. No.                     6.33e+06\n",
            "==============================================================================\n",
            "\n",
            "Notes:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "[2] The condition number is large, 6.33e+06. This might indicate that there are\n",
            "strong multicollinearity or other numerical problems.\n"
          ]
        }
      ],
      "source": [
        "# here\n",
        "# 2 continuous variable with interaction\n",
        "mod = smf.ols(formula='median_house_value ~ total_bedrooms * population\t', data=df)\n",
        "\n",
        "res = mod.fit()\n",
        "\n",
        "print(res.summary())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "652b1c3f",
      "metadata": {
        "id": "652b1c3f"
      },
      "source": [
        "    \n",
        "6. Reading these **multiple linear regression** model specifications just requires understanding straight up simple math... the equations mean exactly what they say and should be interpreted as such; however, that said...\n",
        "\n",
        "    - One thing to be sure to keep in mind when you're interpreting the equations, though, is that you interpret them with respect to one variable at a time, under the assumption \"with all other predictor variables held constant\"\n",
        "    - This \"with all other predictor variables held constant\" assumption is \"theoretical\" in the sense that in observed data you probably can't \"make a one-unit change to just one predictor variable while simultanesouly holding all other predictor variables constant\"\n",
        "    - But, from the perspective of the mathematical equation, you can indeed \"theoretically\" consider this possibility, and this is indeed the perspective that is used to interpret the relationship of each single predictor variable with the outcome on a one-by-one, case-by-case basis\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "863c55e4",
      "metadata": {
        "id": "863c55e4"
      },
      "outputs": [],
      "source": [
        "# ANYWAY, hopefully you're TA picks an interesting data set to demo!\n",
        "# Students *could* potentially help by suggesting an intersting one to try out, too..."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "90f5565d",
      "metadata": {
        "id": "90f5565d"
      },
      "source": [
        "## 💬 🗣️ Communication Activity #2<br>[in the (final 10 minutes of) time remaining<br>(or more or less depending on the time needed for the demo above...)]\n",
        "\n",
        "Return to you <u>**SIX**</u> groups of the first **Communication Activity**, ideally those corresponding to the **Course Project teams** of the TUT\n",
        "\n",
        "1. Discuss **Individual Project Proposals** that might relate to **multiple linear regression** and work to write out a model specification that you could analyze for the **Canadian Social Connection Survey** dataset variables for you **Course Project**\n",
        "2. Explore possible ideas for **outcome variables** (and related relavent **predictor variables**) that you might examine for your **Course Project**\n",
        "3. Consider which of the kinds of **predictor variables** could (and should) be used for your specification, and if you can imagine other **multiple linear regression** model specifications that might take advantage of the variety of types of **predictor variables** that can be leveragecd in the **multiple linear regression** context.\n",
        "\n",
        "> 1. **For your course project**, you're probably working with categorical variables: look into **logistic regression** for now. That is, *if you're understanding* **multiple linear regression** then you can do **logistic regression**.  And next week we'll consider **classification decision trees** which will open up *even further* options for you if you're working with categorical variables.\n",
        "> 2. **Indicator variables** are likely going to be **VERY** important for you as well, regarding the **predictor variables** that you're interested in.\n",
        "> 3. A simpler way to start with considering potential associations between **categorical variables** is to look at **bar plots** of one **categorical variable** across different levels over *another* **categorical variable**. And you can additionally either (a) look into using the **Fisher exact Test**, or (b) think about how you could design **YOUR OWN** hypothesis testing using **simulation** to creating a **sampling distribution under a null hypothesis assumption** in order to produce a **p-value**.\n",
        ">\n",
        ">\n",
        "> 4. <u>**Or as I ALWAYS SAY, \"better yet\", figure how how to create a (perhaps \"double\") _bootstrapped confidence interval_**.</u>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c133a1e5",
      "metadata": {
        "id": "c133a1e5"
      },
      "outputs": [],
      "source": []
    }
  ],
  "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.11.4"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}