{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Linear Regression\n",
    "\n",
    "Student Heigth and Weight example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.stats.power as smp\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Linear regression with Height and weight"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear regression with height and weight data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data\n",
    "x = np.array([168, 161, 167, 179, 184, 166, 198, 187, 191, 179])             # height data\n",
    "y = np.array([65.5, 58.3, 68.1, 85.7, 80.5, 63.4, 102.6, 91.4, 86.7, 78.9])  # weight data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Today we have WEIGHT as our outcome variable (y)\n",
    "\n",
    "and HEIGHT as our explanatory variable (x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ALWAYS start by visualising your data!\n",
    "\n",
    "# we make a scatter plot:\n",
    "plt.scatter(x,y)\n",
    "plt.ylabel(\"weight\")\n",
    "plt.xlabel(\"height\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear relationship?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets put data into a **pandas dataframe** that we call \"student\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "student = pd.DataFrame({'x': x, 'y': y})  \n",
    "print(student)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we do a **linear regression model**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for now we just run the code, we will understand it later\n",
    "\n",
    "# fit the model\n",
    "fitStudents = smf.ols(formula = 'y ~ x', data=student).fit()  # OBS: use the statsmodels.formula.api library (smf)\n",
    "\n",
    "# Get prediction and confidence intervals\n",
    "x_pred = pd.DataFrame({'x': np.arange(160,200, 1)})\n",
    "pred = fitStudents.get_prediction(x_pred).summary_frame(alpha=0.05)\n",
    "\n",
    "plt.scatter(x,y)\n",
    "plt.plot(x_pred, pred[\"mean\"], color=\"red\")\n",
    "plt.ylabel(\"weight\")\n",
    "plt.xlabel(\"height\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red line is the **regression line**.\n",
    "\n",
    "The regression line is a straight line.\n",
    "\n",
    "Is is defined by a **slope** and an **intercept** with the y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# later today we will be able to estimate values for the parameters: intercept and slope\n",
    "\n",
    "# we will also learn to make a \"Regression Results Table\" with Python (using smf library)\n",
    "print(fitStudents.summary(slim=True))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table includes a lot of information and we will understand almost all of it by the end of today :)\n",
    "\n",
    "We will also be able to estimate **standard errors** and **conficence intervals** for the parameters of the model (also in table above)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Later today we will be able to make a **prediction interval** for \"future\" **observations** (with other x-values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot of prediction interval\n",
    "plt.scatter(x,y)\n",
    "plt.plot(x_pred, pred[\"mean\"], color=\"red\")\n",
    "plt.plot(x_pred, pred[\"obs_ci_lower\"], color=\"lightgrey\")\n",
    "plt.plot(x_pred, pred[\"obs_ci_upper\"], color=\"lightgrey\")\n",
    "plt.ylabel(\"weight\")\n",
    "plt.xlabel(\"height\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and later today we will also be able to make a **confidence interval** for the mean at every value on the x-axis (*confidence interval for the line*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add confidence interval for the line\n",
    "plt.scatter(x,y)\n",
    "plt.plot(x_pred, pred[\"mean\"], color=\"red\")\n",
    "plt.plot(x_pred, pred[\"mean_ci_lower\"], color=\"grey\")\n",
    "plt.plot(x_pred, pred[\"mean_ci_upper\"], color=\"grey\")\n",
    "plt.plot(x_pred, pred[\"obs_ci_lower\"], color=\"lightgrey\")\n",
    "plt.plot(x_pred, pred[\"obs_ci_upper\"], color=\"lightgrey\")\n",
    "plt.ylabel(\"weight\")\n",
    "plt.xlabel(\"height\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other things we will talk about:\n",
    "\n",
    "1) How does correlation relate to the linear regression (slope)?\n",
    "\n",
    "2) How much of the total variation in y is explained by the linear relationship?\n",
    "\n",
    "3) Is there a *significant* relationship between x and y?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Estimating parameters $\\beta_0$ and $\\beta_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulated data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we compute som simulated data, so we know the \"real\" beta_0 and beta_1 (and sigma)\n",
    "np.random.seed(23498)\n",
    "\n",
    "beta_0 = 50\n",
    "beta_1 = 200\n",
    "sigma = 90\n",
    "\n",
    "# choose som random x-values:\n",
    "x = stats.uniform.rvs(size = 20, loc=-2, scale = 6)\n",
    "\n",
    "# simulate y-values from statistical model:\n",
    "y = beta_0 + beta_1*x + stats.norm.rvs(size = 20, loc=0, scale = sigma)\n",
    "\n",
    "plt.scatter(x,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is now our \"Toy Data\" :-)\n",
    "\n",
    "Lets see if we can estimate the parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate estimates beta_0_hat and beta_1_hat:\n",
    "\n",
    "Sxx = np.sum((x - x.mean())**2)\n",
    "beta_1_hat = np.sum((x - x.mean())*(y - y.mean())) / Sxx\n",
    "beta_0_hat = y.mean() - beta_1_hat*x.mean()\n",
    "\n",
    "print(beta_0_hat, beta_1_hat)\n",
    "print(beta_0, beta_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot estimated model with data:\n",
    "plt.scatter(x,y)\n",
    "plt.plot([-2,4], [beta_0_hat + beta_1_hat*(-2), beta_0_hat + beta_1_hat*(4)], color=\"red\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "also plot the \"true\" underlying model (this is *only* possible because we have simulated data - IRL one cannot plot the \"true\" model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x,y)\n",
    "plt.plot([-2,4], [beta_0_hat + beta_1_hat*(-2), beta_0_hat + beta_1_hat*(4)], color=\"red\")\n",
    "plt.plot([-2,4], [beta_0 + beta_1*(-2), beta_0 + beta_1*(4)], color=\"black\", linestyle='--')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also estimate the variance of residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = y - beta_0_hat - beta_1_hat * x\n",
    "sigma_hat = np.sqrt(np.sum(res**2)/(20-2))\n",
    "print(sigma_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare this value to the real underlying sigma, which was 90 :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Variation of $\\beta_0$ and $\\beta_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulating **many** samples from the underlying known model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# repeat with new sample:\n",
    "y = beta_0 + beta_1*x + stats.norm.rvs(size = 20, loc=0, scale = sigma) # y-data with new values of residuals\n",
    "\n",
    "# re-calculate parameters\n",
    "Sxx = np.sum((x - x.mean())**2)\n",
    "beta_1_hat = np.sum((x - x.mean())*(y - y.mean())) / Sxx\n",
    "beta_0_hat = y.mean() - beta_1_hat*x.mean()\n",
    "print(beta_0_hat, beta_1_hat)\n",
    "\n",
    "#plot:\n",
    "plt.scatter(x,y)\n",
    "plt.plot([-2,4], [beta_0_hat + beta_1_hat*(-2), beta_0_hat + beta_1_hat*(4)], color=\"red\")\n",
    "plt.plot([-2,4], [beta_0 + beta_1*(-2), beta_0 + beta_1*(4)], color=\"black\", linestyle='--')\n",
    "plt.ylim([-400,1000])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the cell above several times\n",
    "\n",
    "Observe how the estimated model changes due to fluctuations in data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat for many samples and make histograms of estimated parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 1000\n",
    "b0 = np.zeros(k)\n",
    "b1 = np.zeros(k)\n",
    "for i in range(k):\n",
    "   y = beta_0 + beta_1*x + stats.norm.rvs(size = 20, loc=0, scale = sigma)\n",
    "   beta_1_hat = np.sum((x - x.mean())*(y - y.mean())) / Sxx\n",
    "   beta_0_hat = y.mean() - beta_1_hat*x.mean()\n",
    "   b0[i] = beta_0_hat\n",
    "   b1[i] = beta_1_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(1,2,figsize=(8,4))\n",
    "\n",
    "ax0.hist(b0)\n",
    "ax0.axvline(x=beta_0, color=\"black\", linestyle='--')\n",
    "\n",
    "ax1.hist(b1)\n",
    "ax1.axvline(x=beta_1, color=\"black\", linestyle='--')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4: Estimate standard errors of the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recall the height and weigth data:\n",
    "print(student)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 1 ###  Estimate parameters beta_0_hat and bata_1_hat *manually*\n",
    "\n",
    "Sxx = np.sum((student[\"x\"] - student[\"x\"].mean())**2)\n",
    "\n",
    "beta_1_hat = np.sum((student[\"x\"] - student[\"x\"].mean())*(student[\"y\"] - student[\"y\"].mean())) / Sxx\n",
    "\n",
    "beta_0_hat = student[\"y\"].mean() - beta_1_hat*student[\"x\"].mean()\n",
    "\n",
    "print(beta_0_hat, beta_1_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot:\n",
    "plt.scatter(student[\"x\"], student[\"y\"])  # plot the data\n",
    "plt.plot([160, 200], [beta_0_hat + beta_1_hat*160, beta_0_hat + beta_1_hat*200], color=\"red\")  # plot the line (from 160 to 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 2 ###  Estimate standard error for the parameters *manually* \n",
    "\n",
    "# first we add some columns to the dataset:\n",
    "student[\"y_pred\"] = beta_0_hat + beta_1_hat*student[\"x\"]  # y-values ON the regression line\n",
    "student[\"residuals\"] = student[\"y\"] - student[\"y_pred\"]   # residuals\n",
    "print(student)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RSS = np.sum(student[\"residuals\"]**2)\n",
    "sigma_hat = np.sqrt(RSS/(10-2))\n",
    "\n",
    "se_beta_0_hat = sigma_hat*np.sqrt(1/10 + student[\"x\"].mean()**2 / Sxx)\n",
    "\n",
    "se_beta_1_hat = sigma_hat*np.sqrt(1/Sxx)\n",
    "\n",
    "print(\"Intercept: \", beta_0_hat, se_beta_0_hat)\n",
    "print(\"slope:     \", beta_1_hat, se_beta_1_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### automatic ###  Do it all with inbuilt python function smf.ols\n",
    "fitStudents = smf.ols(formula = 'y ~ x', data=student).fit()  # OBS: use the statsmodels.formula.api library (smf)\n",
    "print(fitStudents.summary(slim=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get parameters\n",
    "print(fitStudents.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get parameter standard errors\n",
    "print(fitStudents.bse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get estimate of sigma^2:\n",
    "print(fitStudents.scale)\n",
    "print(sigma_hat**2) #compare to the one we calculated *manualy*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OBS: \"scale\" relates to **variance** og residuals "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get fitted values:\n",
    "print(fitStudents.fittedvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(compare with the \"y_pred\" values that we have computed above)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 5: Hypothesis test for parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitStudents = smf.ols(formula = 'y ~ x', data=student).fit()  # OBS: use the statsmodels.formula.api library (smf)\n",
    "print(fitStudents.summary(slim=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Warning!! The p-values appear to be zero but they are in reality only rounded numbers\n",
    "\n",
    "# print the pvalues seperately to be sure:\n",
    "print(fitStudents.pvalues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also do the test manually (here only for the intercept):\n",
    "\n",
    "# calculate \"t_obs\"\n",
    "t_obs_int = -119.9581/18.897\n",
    "print(t_obs_int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find corresponding p-value:\n",
    "print(2*stats.t.cdf(t_obs_int, df=10-2))  # obs! df = n - 2  (2 is the number of parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 6: Confidence interval for parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitStudents = smf.ols(formula = 'y ~ x', data=student).fit()  \n",
    "print(fitStudents.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "KAHOOT! (x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also calculate CI *manually*\n",
    "\n",
    "# here only for the intercept:\n",
    "\n",
    "print(-119.9581 - stats.t.ppf(0.975, df=10-2)*18.897)\n",
    "print(-119.9581 + stats.t.ppf(0.975, df=10-2)*18.897)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 7: Confidence interval for the line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recall data and model:\n",
    "plt.scatter(student[\"x\"], student[\"y\"])  # plot the data\n",
    "plt.plot([160, 200], [beta_0_hat + beta_1_hat*160, beta_0_hat + beta_1_hat*200], color=\"red\")  # plot the line (from 160 to 200)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets make predictions of (mean) weight for some other height-values. \n",
    "\n",
    "What is the mean weight for persons of height 170cm, 175cm, 195cm or 205 cm?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate this manually (y = beta0_hat + beta1_hat * x)\n",
    "\n",
    "But we can also use Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make predictions (= mean y-value at a given x-value) and corresponding confidence intervals \n",
    "x_new = pd.DataFrame({'x': [170, 175, 195, 205]})\n",
    "\n",
    "pred = fitStudents.get_prediction(x_new).summary_frame(alpha=0.05) \n",
    "\n",
    "print(pred.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The column \"mean\" is the predicted mean weight for each height in x_new\n",
    "\n",
    "We also get a standard error for the mean weight (mean_se) and a conficence interval (for the specified alpha = 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the new predictions of mean weight:\n",
    "\n",
    "plt.scatter(student[\"x\"], student[\"y\"])\n",
    "plt.plot(student[\"x\"], student[\"y_pred\"], color=\"red\")\n",
    "plt.scatter(x_new, pred[\"mean\"], color=\"black\")\n",
    "plt.scatter(x_new, pred[\"mean_ci_lower\"], color=\"black\", marker=\"_\")\n",
    "plt.scatter(x_new, pred[\"mean_ci_upper\"], color=\"black\", marker=\"_\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the confidence interval i larger for more \"extreme\" heights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the confidence interval for large range of values:\n",
    "\n",
    "x_new = pd.DataFrame({'x': np.arange(155, 205,0.1)})\n",
    "pred = fitStudents.get_prediction(x_new).summary_frame(alpha=0.05) \n",
    "\n",
    "plt.scatter(student[\"x\"],student[\"y\"])\n",
    "plt.plot(x_new, pred[\"mean\"], color=\"red\")\n",
    "plt.plot(x_new, pred[\"mean_ci_lower\"], color=\"grey\")\n",
    "plt.plot(x_new, pred[\"mean_ci_upper\"], color=\"grey\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also calculate this interval *manually* (see Method 5.18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 9: Prediction interval for observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same data as above\n",
    "\n",
    "# now we want the prediction interval for individual (new/future) observations\n",
    "\n",
    "print(pred.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the predictions intervals for a wide range of x-values:\n",
    "plt.scatter(student[\"x\"],student[\"y\"])\n",
    "plt.plot(x_new, pred[\"mean\"], color=\"red\")\n",
    "plt.plot(x_new, pred[\"obs_ci_lower\"], color=\"lightgrey\")\n",
    "plt.plot(x_new, pred[\"obs_ci_upper\"], color=\"lightgrey\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the prediction interval and the confidence interval\n",
    "plt.scatter(student[\"x\"],student[\"y\"])\n",
    "plt.plot(x_new, pred[\"mean\"], color=\"red\")\n",
    "plt.plot(x_new, pred[\"obs_ci_lower\"], color=\"lightgrey\")\n",
    "plt.plot(x_new, pred[\"obs_ci_upper\"], color=\"lightgrey\")\n",
    "plt.plot(x_new, pred[\"mean_ci_lower\"], color=\"grey\")\n",
    "plt.plot(x_new, pred[\"mean_ci_upper\"], color=\"grey\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pernille",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}