{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)\n", "\n", "\n", "(reglin_python_notebook)=\n", "# Regressione lineare con Python" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import arviz as az\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import seaborn as sns\n", "from scipy import stats\n", "import statistics\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "from statsmodels.datasets import get_rdataset\n", "import xarray as xr" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Initialize random number generator\n", "RANDOM_SEED = 8927\n", "rng = np.random.default_rng(RANDOM_SEED)\n", "\n", "plt.style.use(\"bmh\")\n", "plt.rcParams[\"figure.figsize\"] = [10, 6]\n", "plt.rcParams[\"figure.dpi\"] = 100\n", "plt.rcParams[\"figure.facecolor\"] = \"white\"\n", "\n", "sns.set_theme(palette=\"colorblind\")\n", "\n", "%load_ext autoreload\n", "%autoreload 2\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'retina'\n", "%config InlineBackend.figure_format = \"svg\"" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Questo tutorial presenta un esercizio di analisi di regressione utilizzando il linguaggio di programmazione Python. In questo esempio analizzeremo i dati di Sahlins sull'economia delle società \"primitive\". {cite:t}`sahlins2013stone` ha studiato la relazione tra la produzione per lavoratore e la proporzione di consumatori rispetto ai lavoratori nelle famiglie utilizzando i dati raccolti nel villaggio di Mazulu, in Africa centrale. La variabile `acres` rappresenta l'area coltivata da ciascuna famiglia, mentre `consumers` indica il rapporto tra il numero di consumatori e il numero di lavoratori in ogni famiglia. I dati sono contenuti nel DataFrame `Sahlins` del pacchetto R `carData` e possono essere importati in questo Notebook utilizzando la funzione `get_rdataset`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
consumersacres
01.001.71
11.081.52
21.151.29
31.153.09
41.202.21
\n", "
" ], "text/plain": [ " consumers acres\n", "0 1.00 1.71\n", "1 1.08 1.52\n", "2 1.15 1.29\n", "3 1.15 3.09\n", "4 1.20 2.21" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = get_rdataset(\"Sahlins\", package=\"carData\", cache=True).data\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
consumersacres
count20.0020.00
mean1.522.16
std0.350.48
min1.001.29
25%1.272.00
50%1.492.19
75%1.652.37
max2.303.09
\n", "
" ], "text/plain": [ " consumers acres\n", "count 20.00 20.00\n", "mean 1.52 2.16\n", "std 0.35 0.48\n", "min 1.00 1.29\n", "25% 1.27 2.00\n", "50% 1.49 2.19\n", "75% 1.65 2.37\n", "max 2.30 3.09" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.describe().round(2)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Creiamo un diagramma a dispersione." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:35:16.124985\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.scatterplot(x=\"consumers\", y=\"acres\", data=df);" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Aggiungiamo la retta dei minimi quadrati." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:35:27.353068\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.regplot(x=\"consumers\", y=\"acres\", ci=None, data=df);\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Troviamo i coefficienti dei minimi quadrati con la funzione `linregress`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Slope: 0.5163200600920598\n", "Intercept: 1.3756445484797928\n" ] } ], "source": [ "slope, intercept, r_value, p_value, std_err = stats.linregress(df[\"consumers\"], df[\"acres\"])\n", "\n", "# Print the coefficients\n", "print(\"Slope:\", slope)\n", "print(\"Intercept:\", intercept)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Per replicare il risultato precedente usando le equazioni dei minimi quadrati, iniziamo trovando la matrice di varianza-covarianza tra `acres` e `consumers`. È importante notare che il valore della varianza, a differenza della covarianza e della correlazione, dipende dai gradi di libertà." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.22766947, 0.06222526],\n", " [0.06222526, 0.12051684]])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.cov(df[\"acres\"], df[\"consumers\"], ddof=1)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "La covarianza tra `x` e `y` è" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.062225263157894735\n" ] } ], "source": [ "s_xy = np.cov(df[\"acres\"], df[\"consumers\"], ddof=1)[0, 1]\n", "print(s_xy)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "La varianza di `x` è" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.12051684210526312\n" ] } ], "source": [ "v_x = np.cov(df[\"acres\"], df[\"consumers\"], ddof=1)[1, 1]\n", "print(v_x)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo la pendenza della retta di regressione, ovvero il coefficiente $b$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.5163200600920598\n" ] } ], "source": [ "b = s_xy / v_x\n", "print(b)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo l'intercetta della retta di regressione, ovvero il coefficiente $a$." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.3756445484797928\n" ] } ], "source": [ "a = np.mean(df[\"acres\"]) - b * np.mean(df[\"consumers\"])\n", "print(a)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Usiamo Matplotlib per generare un diagramma a dispersione a cui aggiungiamo la retta di regressione con i coefficienti calcolati in precedenza. Il risultato è identico a quello trovato da Seaborn." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "x = df[\"consumers\"]\n", "y = df[\"acres\"]\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:37:00.932717\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, y, \"o\")\n", "plt.plot(x, a + b * x);\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Residui\n", "\n", "Consideriamo ora i residui. Iniziamo a calcolare i valori predetti $\\hat{y}$." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 1.891965\n", "1 1.933270\n", "2 1.969413\n", "3 1.969413\n", "4 1.995229\n", "5 2.046861\n", "6 2.083003\n", "7 2.083003\n", "8 2.113982\n", "9 2.129472\n", "10 2.160451\n", "11 2.186267\n", "12 2.227573\n", "13 2.227573\n", "14 2.227573\n", "15 2.232736\n", "16 2.341163\n", "17 2.423774\n", "18 2.434101\n", "19 2.563181\n", "Name: consumers, dtype: float64\n" ] } ], "source": [ "y_hat = a + b * x\n", "print(y_hat)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Rappresentiamo i valori $\\hat{y}$ con dei quadrati rossi nel grafico che riporta la retta di regressione." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:37:12.549430\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, y, \"o\")\n", "plt.plot(x, a + b * x)\n", "plt.plot(x, y_hat, \"s\", color=\"red\");\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo i residui." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 -0.181965\n", "1 -0.413270\n", "2 -0.679413\n", "3 1.120587\n", "4 0.214771\n", "5 0.213139\n", "6 0.316997\n", "7 0.016997\n", "8 -0.153982\n", "9 -0.039472\n", "10 -0.140451\n", "11 -0.876267\n", "12 -0.057573\n", "13 0.052427\n", "14 0.182427\n", "15 -0.002736\n", "16 0.698837\n", "17 -0.363774\n", "18 0.295899\n", "19 -0.203181\n", "dtype: float64\n" ] } ], "source": [ "e = y - y_hat\n", "print(e)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Controlliamo che $\\sum_i(e_i) = 0$." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.1086244689504383e-15" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum(e)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "I residui corrispondono ai segmenti verticali che, nel diagramma a dispersione, collegano i punti al valore predetto sulla retta di regressione." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:37:21.899448\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, y, \"o\")\n", "plt.plot(x, a + b * x)\n", "plt.plot(x, y_hat, \"s\", color=\"red\")\n", "\n", "plt.vlines(x, y, y - e);\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Correlazione e pendenza della retta di regressione\n", "\n", "Se standardizziamo $X$ e $Y$, la pendenza della retta di regressione\n", "\n", "$$\n", "b = \\frac{S_{x,y}}{Var_{x}}\n", "$$\n", "\n", "diventa identica alla correlazione tra le due variabili in quanto, per definizione, la correlazione è la covarianza dei dati standardizzati e, nuovamente per definizione, la varianza di una variabile standardizzata è uguale a 1.\n", "\n", "Inoltre, quando i dati sono standardizzati, $a$ = 0, in quanto \n", "\n", "$$\n", "a = \\bar{y} - b \\cdot \\bar{x}.\n", "$$\n", "\n", "Per ottenere una dimostrazione numerica, creaimo un diagramma a dispersione con le variabili standardizzate, a cui aggiungeremo la retta di regressione con pendenza uguale a $r$ e intercetta 0." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "zx = (x - np.mean(x)) / np.std(x)\n", "zy = (y - np.mean(y)) / np.std(y)\n" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.3756561199682136\n" ] } ], "source": [ "r_xy = np.corrcoef(x, y)[0, 1]\n", "print(r_xy)\n" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-05-08T07:37:30.082398\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.1, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(zx, zy, \"o\")\n", "plt.plot(zx, 0 + r_xy * zx);\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Verifichiamo." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Slope: 0.37565611996821363\n", "Intercept: 2.7789083832427757e-16\n" ] } ], "source": [ "df1 = pd.DataFrame({\"zx\": zx, \"zy\": zy})\n", "\n", "slope, intercept, r_value, p_value, std_err = stats.linregress(df1[\"zx\"], df1[\"zy\"])\n", "\n", "# Print the coefficients\n", "print(\"Slope:\", slope)\n", "print(\"Intercept:\", intercept)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Bontà di adattamento\n", "\n", "### Scomposizione della devianza\n", "\n", "I modelli lineari sono caratterizzati da una scomposizione della devianza totale\n", "\n", "$$\n", "DEV_{\\text{tot}} = \\sum_i{(y_i - \\bar{y})^2}\n", "$$\n", "\n", "in componenti \"spiegate\" e \"non spiegate\".\n", "\n", "Calcoliamo la devianza totale." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.325719999999999\n" ] } ], "source": [ "ss_tot = np.sum((y - np.mean(y)) ** 2)\n", "print(ss_tot)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo la devianza spiegata." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.6104348806456406\n" ] } ], "source": [ "ss_reg = np.sum((y_hat - np.mean(y)) ** 2)\n", "print(ss_reg)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo la devianza non spiegata." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.7152851193543586\n" ] } ], "source": [ "ss_err = np.sum((y - y_hat) ** 2)\n", "print(ss_err)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Verifichiamo." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.325719999999999\n" ] } ], "source": [ "print(ss_reg + ss_err)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Calcoliamo l'indice di determinazione." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.14111752046957288\n" ] } ], "source": [ "r2 = ss_reg / ss_tot\n", "print(r2)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Ovvero, usando una formula equivalente" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.1411175204695727" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1 - ss_err / ss_tot\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Per verificare, eleviamo al quadrato la quantità `r_value` trovata in precedenza. Infatti, nel caso della regressione bivariata, il coefficiente di determinazione è uguale al coefficiente di correlazione di Pearson innalzato al quadrato." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Coefficiente di determinazione: 0.14111752046957288\n" ] } ], "source": [ "print(\"Coefficiente di determinazione:\", r_value**2)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "L'errore standard della regressione è la deviazione standard dei residui calcolata con l'appropriato numero di gradi di libertà. Nel caso della regressione bivariata, i gradi di libertà sono $n-2$, dove $n$ è il numero di osservazioni. Si dice che si \"perdono\" due gradi di libertà a causa del fatto che abbiamo stimato due coefficienti: $a$ e $b$:\n", "\n", "$$\n", "s_e = \\sqrt{\\frac{\\sum_i e^2}{n - 2}}.\n", "$$\n", "\n", "L'errore standard di regressione ci dà un'indicazione approssimativa della media degli errori residui." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.45431787203787166" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(np.sum(e**2) / (len(df[\"consumers\"]) - 2))\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Nel caso attuale, l'errore residuo medio in valore assoluto è" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.31120830458289295" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(np.abs(e))" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "**Esercizio** Si consideri il seguente set di dati: x = 2, 3, 5, 7; y = 11, 12, 14, 13. Disegnare un grafico a dispersione dei dati. Le variabili sembrano correlate? Calcolare la correlazione tra X e Y per verificarlo.\n", "\n", "Definire il modello di regressione con una notazione matematica. Quale segno ci si aspetta per il coefficiente $\\beta$? Spiegare la risposta. Stimare l'intercetta e la pendenza della retta di regressione utilizzando le formule dei minimi quadrati. Il segno del coefficiente $\\beta$ è consistente con la risposta data in precedenza?\n", "\n", "Aggiungere la retta di regressione al grafico a dispersione. Calcolare l'errore standard della regressione. Interpretare in relazione alla rappresentazione grafica dei dati e della retta di regressione.\n", "\n", "Scrivere l'equazione del modello di regressione in termini dei coefficienti stimati con il metodo dei minimi quadrati. Calcolare il valore previsto della variabile Y per un valore di X = 4.\n", "\n", "Standardizzare i dati. Calcolare la pendenza della retta di regressione calcolata per i dati standardizzati. Confrontare il risultato ottenuto con il coefficiente di correlazione.\n", "\n", "Calcolare la devianza totale, la devianza spiegata e la devianza non spiegata. Trovare il coefficiente di determinazione. Interpretare.\n", "\n", "Ripetere i calcoli utilizzando una funzione Python e confrontare i risultati ottenuti in precedenza con quelli riportati dal software." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Inferenza\n", "\n", "Nell'approccio frequentista, l'inferenza viene effettuata mediante l'utilizzo di statistiche $t$ e $p$-value, che vengono comunemente riportati nei risultati di software come Statsmodels. Tali test si basano su teoremi matematici specifici, che non saranno affrontati in questo corso. Al contrario, qui ci concentreremo sull'inferenza bayesiana del modello di regressione, che offre una trattazione più semplice e intuitiva di questo argomento.\n", "\n", "\n", "\n", "## Watermark" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Last updated: Sat May 06 2023\n", "\n", "Python implementation: CPython\n", "Python version : 3.11.3\n", "IPython version : 8.13.2\n", "\n", "scipy : 1.10.1\n", "seaborn : 0.12.2\n", "matplotlib : 3.7.1\n", "arviz : 0.15.1\n", "numpy : 1.23.5\n", "statsmodels: 0.14.0\n", "xarray : 2023.4.2\n", "pandas : 1.5.3\n", "\n" ] } ], "source": [ "%load_ext watermark\n", "%watermark -n -u -v -iv " ] } ], "metadata": { "kernelspec": { "display_name": "pymc", "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" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }