{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IntroStat Week 9\n",
    "\n",
    "Welcome to the 9th lecture in IntroStat\n",
    "\n",
    "During the lectures we will present both slides and notebooks. \n",
    "\n",
    "This is the notebook used in the lecture in week 9.\n"
   ]
  },
  {
   "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": [
    "### Example: Ozon concentration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Today we will be using the same data for all examples\n",
    "\n",
    "The **Airquality dataset**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the \"airquality\" dataset and print first 10 rows:\n",
    "Air = sm.datasets.get_rdataset(\"airquality\").data\n",
    "print(Air.head(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Daily readings of the following air quality values for May 1, 1973 (a Tuesday) to September 30, 1973:\n",
    "\n",
    "- Ozone: Mean ozone in parts per billion from 1300 to 1500 hours at Roosevelt Island\n",
    "\n",
    "- Solar.R: Solar radiation in Langleys in the frequency band 4000–7700 Angstroms from 0800 to 1200 hours at Central Park\n",
    "\n",
    "- Wind: Average wind speed in miles per hour at 0700 and 1000 hours at LaGuardia Airport\n",
    "\n",
    "- Temp: Maximum daily temperature in degrees Fahrenheit at La Guardia Airport.\n",
    "\n",
    "Notice the dataset has some missing values (NaN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many rows? (obs: number of observations is smaller due to missing values in dataset)\n",
    "print(len(Air))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The dataset is now a *Pandas Dataframe*\n",
    "print(type(Air))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When working with \"real\" data one typically need to do some *data-wrangling* \n",
    "\n",
    "This may include joining data from different tables, re-ordering data, re-structuring the dataset\n",
    "\n",
    "Today we need to rename a column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rename \"Solar.R\", because the use of \".\" means something else in python and this could give rise to errors:\n",
    "Air = Air.rename(columns={\"Solar.R\": \"SolarR\"})\n",
    "print(Air.head(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We are interested in the Ozone values \n",
    "\n",
    "How does Ozone depend on solar radiation, wind and temperature?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise the data\n",
    "\n",
    "# split to make 3 plots:\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16,4))\n",
    "\n",
    "# scatter plots with Ozone versus the three other columns\n",
    "Air.plot.scatter('Wind',    'Ozone', ax=ax0)\n",
    "Air.plot.scatter('SolarR',  'Ozone', ax=ax1)\n",
    "Air.plot.scatter('Temp',    'Ozone', ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second plot seems to have variance increasing with increasing solar radiation\n",
    "\n",
    "The dataset seems to have more Ozone values that are relatively small and fewer values that are relatively large\n",
    "\n",
    "We will now inspect the distribution of Ozone values (and consider whether a data-transformation a appropriate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do a histogram of \"Ozone\" column\n",
    "plt.hist(Air[\"Ozone\"], density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since the distribution is very skewed we will do a data-transformation and compute the logarithm of the Ozone column:\n",
    "Air[\"logOzone\"] = np.log(Air[\"Ozone\"])   # this creates an axtra column in the dataframe\n",
    "print(Air)\n",
    "\n",
    "plt.hist(Air[\"logOzone\"], density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This distribution is nicer (but maybe not perfect - the value with logOzone=0 could be considered an outlier)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise the transformed data\n",
    "\n",
    "# split to make 3 plots:\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16,4))\n",
    "\n",
    "# scatter plots with logOzone versus the three other columns\n",
    "Air.plot.scatter('Wind',    'logOzone', ax=ax0)\n",
    "Air.plot.scatter('SolarR',  'logOzone', ax=ax1)\n",
    "Air.plot.scatter('Temp',    'logOzone', ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considerations:\n",
    "\n",
    "The datapoint with logOzone=0 could be an outlier - maybe we should remove it from the dataset.\n",
    "\n",
    "We could model logOzone as a linear function of either wind, solar radiation or temperature. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the outlier:\n",
    "Air = Air[(Air[\"logOzone\"] != 0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(Air))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After removing outlier we have **152 rows** in the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now we do **Simpel Linear Regression** with one explanatory variable at the time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# logOzone and temperature:\n",
    "fit_temp = smf.ols(formula = 'logOzone ~ Temp', data=Air).fit()\n",
    "print(fit_temp.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **>>> Kahoot 1 (3 questions)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get predictions based on the model:\n",
    "temp_range = pd.DataFrame({'Temp': np.arange(50,110, 10)})\n",
    "print(temp_range)\n",
    "\n",
    "pred = fit_temp.get_prediction(temp_range).summary_frame(alpha=0.05)\n",
    "print(pred.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **>>> Kahoot 2 (1 questions)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot data with model predictions:\n",
    "Air.plot.scatter('Temp',    'logOzone')\n",
    "plt.plot(temp_range, pred[\"mean\"], color=\"red\")\n",
    "plt.scatter(temp_range, pred[\"mean\"], color=\"red\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice - this is a *simple linear regression model* with temperature as explanatory variable (and logOzone as dependent variable / outcome variable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets try another explanatory variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# logOzone verus solar radiation:\n",
    "\n",
    "# fit the model\n",
    "fit_rad = smf.ols(formula = 'logOzone ~ SolarR', data=Air).fit()\n",
    "print(fit_rad.summary(slim=True))\n",
    "\n",
    "# Get prediction based on the model:\n",
    "rad_range = pd.DataFrame({'SolarR': np.arange(0,360, 50)})\n",
    "pred = fit_rad.get_prediction(rad_range).summary_frame(alpha=0.05)\n",
    "\n",
    "# Plot data with model:\n",
    "Air.plot.scatter('SolarR',    'logOzone')\n",
    "plt.plot(rad_range, pred[\"mean\"], color=\"red\")\n",
    "plt.scatter(rad_range, pred[\"mean\"], color=\"red\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good (?)\n",
    "\n",
    "Now we try the last explanatory variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# logOzone verus wind:\n",
    "\n",
    "# fit the model\n",
    "fit_wind = smf.ols(formula = 'logOzone ~ Wind', data=Air).fit()\n",
    "print(fit_wind.summary(slim=True))\n",
    "\n",
    "# Get prediction based on the model:\n",
    "wind_range = pd.DataFrame({'Wind': np.arange(0,30,5)})\n",
    "pred = fit_wind.get_prediction(wind_range).summary_frame(alpha=0.05)\n",
    "\n",
    "# Plot data with model predictions:\n",
    "Air.plot.scatter('Wind',    'logOzone')\n",
    "plt.plot(wind_range, pred[\"mean\"], color=\"red\")\n",
    "plt.scatter(wind_range, pred[\"mean\"], color=\"red\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Can we somehow combine all the variables into one model?? \n",
    "\n",
    "#### The answer is yes - with **\"Multiple Linear Regression\"** \n",
    "\n",
    "(return to slides/blackboard before proceding)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Ozon concentration with Multiple Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mulitiple Linear Regression\n",
    "\n",
    "# fit the full model with 3 explanatopry variables:\n",
    "fit_full_model = smf.ols(formula = 'logOzone ~ Temp + SolarR + Wind', data=Air).fit()  # Notice how the model is formulated in python!\n",
    "print(fit_full_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a coefficient (parameter estimate) for each explanatory variable in the model (+ for the intercept)\n",
    "\n",
    "Each parameter estimate also has a Standard Error\n",
    "\n",
    "For each parameter a t-test is performed (t_obs and p-values and confidence intervals are calculated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Confidence and prediction intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_new    = pd.DataFrame({'Temp':   [70],   # a bit of \"wrangling\" is needed\n",
    "                         'SolarR': [200],\n",
    "                         'Wind':   [10]})\n",
    "print(x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = fit_full_model.get_prediction(x_new).summary_frame(alpha=0.05)    # get predictions\n",
    "print(pred.head(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have predicted outcome for given values of explanatory variables \n",
    "\n",
    "(Make sure you understand the meaning of all the printed values above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will now predict for MANY values of Temp and Wind (but keep SolarR constant in order to be able to plot in 3D)\n",
    "\n",
    "Temp_range = np.linspace(50, 100, 100)\n",
    "Wind_range = np.linspace(0, 25, 100)\n",
    "SolarR_new = Air[\"SolarR\"].mean()\n",
    "\n",
    "Temp_new, Wind_new = np.meshgrid(Temp_range, Wind_range)        # make a meshgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "New_data = pd.DataFrame({'Temp':   Temp_new.reshape(1,-1)[0],   # a bit of \"wrangling\" is needed\n",
    "                         'SolarR': SolarR_new,\n",
    "                         'Wind':   Wind_new.reshape(1,-1)[0]})\n",
    "print(New_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred = fit_full_model.get_prediction(New_data).summary_frame(alpha=0.05)    # get predictions\n",
    "print(pred.head(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "New_data_with_pred = pd.concat([New_data, pred], axis=1)                    # merge \"x\"-data with predictions\n",
    "print(New_data_with_pred.head(5))\n",
    "print(len(New_data_with_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 3D plot \n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# Plot a surface\n",
    "ax.plot_surface(np.array(New_data_with_pred[\"Temp\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"Wind\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"mean\"]).reshape(100,100))\n",
    "\n",
    "ax.plot_surface(np.array(New_data_with_pred[\"Temp\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"Wind\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"mean_ci_upper\"]).reshape(100,100), color=\"grey\")\n",
    "\n",
    "ax.plot_surface(np.array(New_data_with_pred[\"Temp\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"Wind\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"mean_ci_lower\"]).reshape(100,100), color=\"grey\")\n",
    "\n",
    "ax.plot_surface(np.array(New_data_with_pred[\"Temp\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"Wind\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"obs_ci_upper\"]).reshape(100,100), color=\"lightgrey\")\n",
    "\n",
    "ax.plot_surface(np.array(New_data_with_pred[\"Temp\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"Wind\"]).reshape(100,100),\n",
    "                np.array(New_data_with_pred[\"obs_ci_lower\"]).reshape(100,100), color=\"lightgrey\")\n",
    "\n",
    "# Add labels\n",
    "ax.set_xlabel('Temp')\n",
    "ax.set_ylabel('Wind')\n",
    "ax.set_zlabel('logOzone')\n",
    "\n",
    "ax.view_init(elev=30, azim=20)\n",
    "#ax.view_init(elev=20, azim=10)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The blue plane are the predicted values. \n",
    "\n",
    "The light grey is the prediction interval.\n",
    "\n",
    "The dark grey is the confidence-interval (for the plane)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(back to slides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Model control\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# full model with 3 explanatory variables:\n",
    "fit_full_model = smf.ols(formula = 'logOzone ~ Temp + SolarR + Wind', data=Air).fit()\n",
    "\n",
    "# residuals:\n",
    "residuals = fit_full_model.resid\n",
    "fittedvalues = fit_full_model.fittedvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot residuals\n",
    "fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,4))\n",
    "\n",
    "# qq-plot of resiudals:\n",
    "sm.qqplot(residuals,ax=ax1, line='s')\n",
    "ax1.set_title(\"QQ plot of residuals\")\n",
    "\n",
    "# plot residuals versus fitted values:\n",
    "ax2.scatter(fittedvalues,residuals)\n",
    "ax2.set_xlabel(\"fitted values\")\n",
    "ax2.set_ylabel(\"residuals\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "qq-plot: check normality of residuals (this one looks fine because the points are almost on the straight line)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scatter plot: is used to look for patterns in residuals (here we se no obvious pattern - at least not a strong one)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot residuals versus individual explanatory variables:\n",
    "Air[\"residuals\"] = residuals\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1,3,figsize=(12,4))\n",
    "\n",
    "# plot residuals versus fitted values:\n",
    "ax0.scatter(Air[\"Temp\"],Air[\"residuals\"])\n",
    "ax0.set_xlabel(\"Temp\")\n",
    "ax0.set_ylabel(\"residuals\")\n",
    "\n",
    "ax1.scatter(Air[\"SolarR\"],Air[\"residuals\"])\n",
    "ax1.set_xlabel(\"SolarR\")\n",
    "ax1.set_ylabel(\"residuals\")\n",
    "\n",
    "ax2.scatter(Air[\"Wind\"],Air[\"residuals\"])\n",
    "ax2.set_xlabel(\"Wind\")\n",
    "ax2.set_ylabel(\"residuals\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see no obvious patterns in the residuals - which is good!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: curvilinear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the original data (but with log(ozone))\n",
    "\n",
    "# split to make 4 plots:\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16,4))\n",
    "\n",
    "# scatter plots with logOzone versus the three other columns\n",
    "Air.plot.scatter('Wind',    'logOzone', ax=ax0)\n",
    "Air.plot.scatter('SolarR', 'logOzone', ax=ax1)\n",
    "Air.plot.scatter('Temp',    'logOzone', ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Would it be a good idea to include some squared terms?\n",
    "\n",
    "Lets try with **Wind**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linear model with wind (simple linear regression):\n",
    "\n",
    "# fit the model\n",
    "fit_wind = smf.ols(formula = 'logOzone ~ Wind', data=Air).fit()\n",
    "print(fit_wind.summary(slim=True))\n",
    "\n",
    "# Get prediction based on the model:\n",
    "wind_range = pd.DataFrame({'Wind': np.arange(0,25.5,.5)})\n",
    "pred = fit_wind.get_prediction(wind_range).summary_frame(alpha=0.05)\n",
    "\n",
    "# Plot data with model predictions:\n",
    "Air.plot.scatter('Wind',    'logOzone')\n",
    "plt.plot(wind_range, pred[\"mean\"], color=\"red\")\n",
    "plt.scatter(wind_range, pred[\"mean\"], color=\"red\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make column with wind-squared /values:\n",
    "Air[\"Wind2\"] = Air[\"Wind\"]**2\n",
    "\n",
    "print(Air.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model with Wind and Wind2\n",
    "\n",
    "# fit the model\n",
    "fit_wind2 = smf.ols(formula = 'logOzone ~ Wind + Wind2', data=Air).fit()\n",
    "print(fit_wind2.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks great - wind-squared is also significant. \n",
    "\n",
    "Is it a better model than the linear one?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visual inspection:\n",
    "\n",
    "# Get prediction based on the model:\n",
    "wind_range[\"Wind2\"] = wind_range[\"Wind\"]**2\n",
    "pred2 = fit_wind2.get_prediction(wind_range).summary_frame(alpha=0.05)\n",
    "\n",
    "# Plot data with model predictions:\n",
    "Air.plot.scatter('Wind',    'logOzone')\n",
    "plt.plot(wind_range[\"Wind\"], pred2[\"mean\"], color=\"red\")\n",
    "plt.scatter(wind_range[\"Wind\"], pred2[\"mean\"], color=\"red\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is it a better model?\n",
    "\n",
    "Is it an *\"overfit\"*?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets have a quick look at the residuals (model control) in both linear and quadratic cases:\n",
    "# residuals:\n",
    "residuals = fit_wind.resid\n",
    "fittedvalues = fit_wind.fittedvalues\n",
    "\n",
    "residuals2 = fit_wind2.resid\n",
    "fittedvalues2 = fit_wind2.fittedvalues\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(2,3,figsize=(14,8))\n",
    "\n",
    "# for linear model:\n",
    "# qq-plot of resiudals:\n",
    "sm.qqplot(residuals,ax=axs[0,0], line='s')\n",
    "axs[0,0].set_title(\"QQ plot of residuals - Linear Model\")\n",
    "\n",
    "# plot residuals versus fitted values:\n",
    "axs[0,1].scatter(fittedvalues,residuals)\n",
    "axs[0,1].set_xlabel(\"fitted values\")\n",
    "axs[0,1].set_ylabel(\"residuals\")\n",
    "axs[0,1].set_title(\"Residuals - Linear Model\")\n",
    "\n",
    "# plot data with model\n",
    "Air.plot.scatter('Wind',    'logOzone', ax = axs[0,2])\n",
    "axs[0,2].plot(wind_range[\"Wind\"], pred[\"mean\"], color=\"red\")\n",
    "\n",
    "# for quadratic model:\n",
    "# qq-plot of resiudals:\n",
    "sm.qqplot(residuals2,ax=axs[1,0], line='s')\n",
    "axs[1,0].set_title(\"QQ plot of residuals - Quadratic Model\")\n",
    "\n",
    "# plot residuals versus fitted values:\n",
    "axs[1,1].scatter(fittedvalues2,residuals2)\n",
    "axs[1,1].set_xlabel(\"fitted values\")\n",
    "axs[1,1].set_ylabel(\"residuals\")\n",
    "axs[1,1].set_title(\"Residuals - Quadratic Model\")\n",
    "\n",
    "# plot data with model\n",
    "Air.plot.scatter('Wind',    'logOzone', ax = axs[1,2])\n",
    "axs[1,2].plot(wind_range[\"Wind\"], pred2[\"mean\"], color=\"red\")\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which model is better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **>>> Kahoot (1 question)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe it would make more sense to include a quadratic term for the SolarR-variable. You can try this yourself! (simply change Wind and Wind2 to SolarR and SolarR2 in the code above)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Colinearity and confounding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets make a FAKE variable which is super CORRELATED TO WIND:\n",
    "Air[\"fake_wind\"] = 2*Air[\"Wind\"] + stats.norm.rvs(size=len(Air), loc=0, scale=1)\n",
    "print(Air.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Air.plot.scatter('fake_wind',    'logOzone')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Air.plot.scatter('fake_wind',    'Wind')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we fit a model using both Wind and fake_wind:\n",
    "fit_full_model = smf.ols(formula = 'logOzone ~ Temp + SolarR + Wind + fake_wind', data=Air).fit()\n",
    "print(fit_full_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both Wind and fake_wind become insignificant!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **>>> Kahoot (1 question)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Model selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# re-visit the Multiple Linear Regression model:\n",
    "print(fit_full_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we include even more explanatory variables?\n",
    "\n",
    "Are more variables always better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How do we select the best model? How do we choose which explanatory variables to include?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First imagine we had even MORE explanatory variables:\n",
    "\n",
    "# lets make a \"fake\" variable (that is independent of the other variables):\n",
    "Air[\"fake\"] = stats.norm.rvs(size=len(Air))\n",
    "print(Air.head())\n",
    "\n",
    "Air.plot.scatter('fake',    'logOzone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would happen if we added this \"fake\" variable to the full model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_full_model = smf.ols(formula = 'logOzone ~ Temp + SolarR + Wind + fake', data=Air).fit()\n",
    "print(fit_full_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"slope\" parameter for \"fake\" is not significant (not significantly different from zero). We see this from the p-value, which is larger than 0.05.\n",
    "\n",
    "This is because \"fake\" *carries no information* - it is just random values and has no correlation with the outcome variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Model Selection:**\n",
    "\n",
    "So the question is how to choose which variable to include in a model - its not always better to include more variables. \n",
    "\n",
    "There are different approaches for selecting which explanatory variables to include in the model (*model selection*):\n",
    "\n",
    "*Backward selection*: Remove the least significant variables one at the time\n",
    "\n",
    "*Forward selection*: Start with the most significant variable and add one at the time\n",
    "\n",
    "OBS: There is no perfect way to do model selection! One should also be aware of possible correlations in explanatory variables (collinearity). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: 3 treatment groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treatment_data = pd.read_csv(\"treatment_data.csv\", sep=',')  # for this code to work, you need to save the treatment_data.csv file in the same folder as you keep this notebook\n",
    "print(treatment_data)\n",
    "print(len(treatment_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable \"Treat\" is qualitative and indicates which treatment group the patient belongs to. \n",
    "\n",
    "The variables Prewt and Postwt are is the weight of the patient before and after treatment. \n",
    "\n",
    "The variable Diff is simply (Postwt - Prewt)\n",
    "\n",
    "OBS: The treatment is for patients with eating disorder - the goal of the treatment is to make patients gain weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many values in the categorical column (the column with text)?\n",
    "print(treatment_data[\"Treat\"].value_counts())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise data (here we group by treatment type)\n",
    "\n",
    "treatment_data.boxplot(\"diff\", by='Treat', showmeans=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean values (of each group) are indicated by triangles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare two groups with **linear regression**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sample data with only two of the treatment groups:\n",
    "sample_data = treatment_data[treatment_data['Treat'].isin([\"Cont\", \"FT\"])]\n",
    "print(len(sample_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_data.boxplot(\"diff\", by='Treat', showmeans=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could of course compare these two groups with the methods we have learned for 2 samples. \n",
    "\n",
    "The comparison can be formulated as a linear regression problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a model with categorical variable:\n",
    "sample_model = smf.ols(formula = 'diff ~ Treat', data=sample_data).fit()\n",
    "print(sample_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the coefficient for FT treatment is 3.4962, which corresponds to the difference between the two means (see boxplot above).\n",
    "\n",
    "Here the \"Cont\" group is chosen as a \"reference group\" and the effect of \"FT\" is shown in the model (Treat[T.FT] = 3.4962). This means that the average \"diff\" (i.e. the average weight-gain) in the FT group is 3.4962 higher than in the Cont group. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple linear regression with categorical AND quantitative variables\n",
    "\n",
    "We can also add Prewt in the regression model - maybe the weight-gain (\"diff\") also depends on the weight prior to treatment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot diff versus Prewt (pre-treatment weight)\n",
    "# and indicate treatment group by color\n",
    "colors = {'Cont':'black', 'CBT':'red', 'FT':'blue'}\n",
    "treatment_data.plot.scatter(\"Prewt\", \"diff\", c=treatment_data[\"Treat\"].map(colors))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that diff (weight-gain) is higher for the persons that had the lowest starting weight (smallest Prewt values).\n",
    "\n",
    "The colors indicate the different treatment groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model with categorical variabel (+ continous variabel):\n",
    "cat_model = smf.ols(formula = 'diff ~ Prewt + Treat', data=treatment_data).fit()\n",
    "print(cat_model.summary(slim=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "x_values = np.arange(32,47,1)\n",
    "y_values_cont = 22.5474 - 0.5655 * x_values - 1.8575\n",
    "y_values_cbt  = 22.5474 - 0.5655 * x_values + 0\n",
    "y_values_ft   = 22.5474 - 0.5655 * x_values + 2.0665\n",
    "\n",
    "colors = {'Cont':'black', 'CBT':'red', 'FT':'blue'}\n",
    "treatment_data.plot.scatter(\"Prewt\", \"diff\", c=treatment_data[\"Treat\"].map(colors))\n",
    "plt.plot(x_values, y_values_cont, color=\"black\")\n",
    "plt.plot(x_values, y_values_cbt, color=\"red\")\n",
    "plt.plot(x_values, y_values_ft, color=\"blue\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a regression model with \"Prewt\" and \"treatment group\" as explanatory variables. \n",
    "\n",
    "The effect of Prewt is the same in all groups (the slope is the same in all groups, by design) \n",
    "\n",
    "The 3 groups have different \"levels\" "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Pernille_Env_python_311",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}