D. Malchiodi, Superhero data science. Vol 1: probabilità e statistica: Analisi della varianza.


Analisi della varianza

Ipotizziamo di avere a disposizione delle osservazioni di un medesimo attributo divise in $G$ gruppi, per esempio perché si tratta del reddito di individui che vivono in diverse città, oppure di un valore ematico di pazienti sottoposti a diversi trattamenti clinici e così via. Formalmente, indichiamo rispettivamente con $n_1, \dots, n_G$ le numerosità dei vari gruppi, con $n = n_1 + \dots + n_G$ il numero totale di osservazioni e, fissato $g \in \{1, \dots, G\}$ e $i \in \{1, \dots, n_g\}$, denotiamo con $x^g_i$ il valore dell'$i$-esima osservazione nel gruppo $g$.

Se si è interessati a valutare l'ipotesi che i valori delle medie nei vari gruppi siano sensibilmente differenti, per esempio perché si vuole dimostrare che il reddito non sia troppo diverso in un gruppo di città, oppure per dimostrare l'efficacia di un dato trattamento medico, è possibile applicare un metodo chiamato ANOVA (ANalysis Of VAriance). L'idea alla base di questo metodo è che se non vi sono sostanziali differenze tra i gruppi considerati, allora calcolare la varianza all'interno di un gruppo qualsiasi non dovrebbe portare a un risultato molto dissimile da quello ottenuto effettuando il calcolo su tutti i dati a disposizione. Più formalmente, definite:

  • la media campionaria $\overline x = \frac{1}{n} \sum_{g=1}^G \sum_{i=1}^{n_g} x^g_i$ su tutte le osservazioni;

  • la media campionaria $\overline x^g = \frac{1}{n_g} \sum_{i=1}^{n_g} x^g_i$ all'interno del $g$-esimo gruppo, per ogni $g = 1, \dots, n_G$;

  • la somma totale degli scarti

\begin{equation} \mathrm{SS}_{\mathrm T} = \sum_{g=1}^G \sum_{i=1}^{n_g} \left( x^g_i - \overline x \right)^2; \end{equation}
  • la somma degli scarti entro i gruppi (o, usando la terminologia inglese, within groups)
\begin{equation} \mathrm{SS}_{\mathrm W} = \sum_{g=1}^G \sum_{i=1}^{n_g} \left( x^g_i - \overline x^g \right)^2; \end{equation}
  • la somma degli scarti tra i gruppi (o, usando la terminologia inglese, between groups), pesata rispetto alla numerosità dei vari gruppi:
\begin{equation} \mathrm{SS}_{\mathrm B} = \sum_{g=1}^G n_g \left( \overline x^g - \overline x \right)^2; \end{equation}

A partire da ognuna di queste somme è facile calcolare le corrispondenti varianze campionarie:

  • la varianza campionaria su tutte le osservazioni: $s^2_{\mathrm T} = \frac{1}{n-1} \mathrm{SS}_{\mathrm T}$;

  • la varianza campionaria delle medie tra i gruppi: $s^2_{\mathrm B} = \frac{1}{G-1} \mathrm{SS}_{\mathrm B}$ (il motivo per cui viene fatta la divisione per $G-1$ è analogo alla ragione per cui il calcolo della varianza campionaria viene fatta dividendo per $n-1$, e richiede un maggiore approfondimento teorico per poter essere giustificato);

  • la varianza campionaria dei valori entro i gruppi: $s^2_{\mathrm W} = \frac{1}{n-G} \mathrm{SS}_{\mathrm W}$.

Si può mostrare (chi è interessato può leggere il paragrafo opzionale che segue) che $\mathrm{SS}_{\mathrm T} = \mathrm{SS}_{\mathrm W} + \mathrm{SS}_{\mathrm B}$, e quindi che

\begin{align} \frac{\mathrm{SS}_{\mathrm T}}{n-1} &= \frac{\mathrm{SS}_{\mathrm W}}{n-1} +\frac{\mathrm{SS}_{\mathrm B}}{n-1}, \\ &= \frac{n-G}{n-1} \frac{\mathrm{SS}_{\mathrm W}}{n-G} + \frac{G-1}{n-1} \frac{\mathrm{SS}_{\mathrm B}}{G-1}. \end{align}

Pertanto,

\begin{equation} \text{varianza totale} = \frac{n-G}{n-1} \text{varianza entro i gruppi} + \frac{G-1}{n-1} \text{varianza tra i gruppi}. \end{equation}

Possiamo utilizzare questa uguaglianza per validare o confutare l'ipotesi che le medie nei gruppi siano diverse: se la varianza totale e la varianza entro i gruppi assumono valori non troppo diversi tra loro (e dunque la varianza tra i gruppi risulta trascurabile), allora si confuta l'ipotesi; se al contrario la varianza tra i gruppi assume un valore elevato si può convalidare l'ipotesi.

In realtà per validare o confutare questi tipi di ipotesi esistono metodi quantitativi basati sul calcolo di indici che coinvolgono le quantità che noi abbiamo confrontato in modo qualitativo. L'uso di questi metodi richiede però un approfondimento teorico che va al di là del programma del nostro corso.

Per applicare la tecnica di analisi della varianza selezionaniamo per esempio i due gruppi dei supereroi corrispondenti alle due case editrici più rappresentate: Marvel e DC, e concentriamoci per ogni gruppo sull'indice di forza.

In [1]:
%matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.constants import golden

plt.style.use('fivethirtyeight')
plt.rc('figure', figsize=(5.0, 5.0/golden))


heroes = pd.DataFrame.from_csv('data/heroes.csv', sep=';')

marvel_strength = heroes[(heroes['Publisher'] == 'Marvel Comics') &
                         (pd.notnull(heroes['Strength']))]['Strength']
dc_strength = heroes[(heroes['Publisher'] == 'DC Comics') & \
                     (pd.notnull(heroes['Strength']))]['Strength']
La cella precedente mostra come effettuare un'operazione di selezione complessa su di un _dataframe_, basandosi sull'operatore `&` di congiungzione logica. Analogamente, l'operatore `|` permette di calcolare disgiunzioni logiche.

Iniziamo calcolando $\mathrm{SS}_{\mathrm T}$.

In [2]:
all_strength = pd.concat([marvel_strength, dc_strength])
sum_total = sum((all_strength - all_strength.mean())**2)
sum_total
Out[2]:
561319.43897637713

Analogamente calcoliamo $\mathrm{SS}_{\mathrm W}$.

In [3]:
sum_within = sum((marvel_strength - marvel_strength.mean())**2) + \
                    sum((dc_strength - dc_strength.mean())**2)
sum_within
Out[3]:
561127.62247600162

Infine, calcoliamo $\mathrm{SS}_{\mathrm B}$.

In [4]:
sum_between = len(marvel_strength) * (marvel_strength.mean() - all_strength.mean())**2 + \
                    len(dc_strength) * (dc_strength.mean() - all_strength.mean())**2
sum_between
Out[4]:
191.81650037662837

Verifichiamo innanzitutto che valga, nei limiti dell'approssimazione in virgola mobile, l'uguaglianza $\mathrm{SS}_{\mathrm T} = \mathrm{SS}_{\mathrm W} + \mathrm{SS}_{\mathrm B}$.

In [5]:
sum_total - sum_within - sum_between
Out[5]:
-1.1153815648867749e-09

Calcoliamo infine la varianza totale e la varianza entro i gruppi utilizzando le formule sopra descritte.

In [6]:
n = len(all_strength)

total_var = sum_total / float(n-1)
total_var
Out[6]:
1107.1389328922626
In [7]:
within_var = sum_within / float(n-2)
within_var * (n-2) / (n-1)
Out[7]:
1106.7605965996088

I due valori sono molto vicini, e quindi possiamo avvalorare l'ipotesi che i valori medi dell'indice di forza siano sostanzialmente uguali.

Riscriviamo il codice in modo da incorporare in una funzione la procedura di analisi della varianza, in modo da gestire anche più di due gruppi di osservazioni: la funzione accetterà una lista di tali gruppi come argomento, e restituirà una coppia i cui elementi saranno rispettivamente la varianza totale e la varianza entro i gruppi.

In [8]:
def anova(groups):
    all_elements = pd.concat(groups)
    
    sum_total = sum((all_elements - all_elements.mean())**2)
    sum_within = sum([sum((g - g.mean())**2) for g in groups])
    
    sum_between = sum([len(g) * (g.mean()-all_elements.mean())**2 for g in groups])
    assert(np.abs(sum_total - sum_within - sum_between) < 10**-5)
    n = len(all_elements)
    total_var = sum_total / float(n-1)
    within_var = sum_within / float(n-len(groups))
    
    return (total_var, within_var*(n-len(groups))/(n-1))

Verifichiamo che i valori restituti per i due gruppi già considerati siano gli stessi.

In [9]:
anova([dc_strength, marvel_strength])
Out[9]:
(1107.1389328922639, 1106.7605965996088)

Applicando la procedura all'indice di forza suddiviso tra supereroi e supereroine si ottengono due valori tutto sommato relativamente simili, così che non si possa avvalorare l'ipotesi che l'indice di forza sia distribuito in modo sostanzialmente diverso tra i due generi.

In [10]:
male_strength = heroes[(heroes['Gender'] == 'M') & (pd.notnull(heroes['Strength']))]['Strength']
female_strength = heroes[(heroes['Gender'] == 'F') & (pd.notnull(heroes['Strength']))]['Strength']
anova([male_strength, female_strength])
Out[10]:
(1083.4610943430719, 1072.1108493522443)

Le cose cambiano se consideriamo la divisione tra generi per i supereroi DC, valutando la differenza nella distribuzione del peso.

In [11]:
male_year = heroes[(heroes['Publisher'] == 'DC Comics') & \
                   (heroes['Gender'] == 'M') & \
                   (pd.notnull(heroes['Weight']))]['Weight']
female_year = heroes[(heroes['Publisher'] == 'DC Comics') & \
                     (heroes['Gender'] == 'F') & \
                     (pd.notnull(heroes['Weight']))]['Weight']
anova([male_year, female_year])
Out[11]:
(9314.5881418390782, 9049.0126309058724)

Dimostrazione *

Innanzitutto notiamo che per ogni $g = 1, \dots, G$,

\begin{equation} \sum_{i=1}^{n_g} x^g_i = n_g \overline x^g, \end{equation}

e quindi

\begin{align} \mathrm{SS}_{\mathrm T} &= \sum_{g=1}^G \sum_{i=1}^{n_g}\left( \left(x^g_i\right)^2 - 2 \overline x x^g_i + \left( \overline x \right)^2 \right) \\ &= \sum_{g=1}^G \sum_{i=1}^{n_g}\left( \left(x^g_i\right)^2 - 2 \overline x x^g_i + \left( \overline x \right)^2 + \left( \overline x^g \right)^2 - \left( \overline x^g \right)^2 + 2 x^g_i \overline x^g - 2 x^g_i \overline x^g \right) = \\ &= \sum_{g=1}^G \sum_{i=1}^{n_g} \left( x^g_i - \overline x^g \right)^2 + \sum_{g=1}^G \sum_{i=1}^{n_g} \left( \left( \overline x \right)^2 - \left( \overline x^g \right)^2 -2 \overline x x^g_i + 2 x^g_i \overline x^g \right) = \\ &= \mathrm{SS}_{\mathrm W} + \sum_{g=1}^G \left( n_g \left( \overline x \right)^2 - n_g \left( \overline x^g \right)^2 - 2 \overline x \sum_{i=1}^{n_g} x^g_i + 2 \overline x^g \sum_{i=1}^{n_g} x^g_i \right) = \\ &= \mathrm{SS}_{\mathrm W} + \sum_{g=1}^G n_g \left( \left( \overline x \right)^2 - \left( \overline x^g \right)^2 - 2 \overline x \overline x^g + 2 \left( \overline x^g \right)^2 \right) = \\ &= \mathrm{SS}_{\mathrm W} + \sum_{g=1}^G n_g \left( \overline x^g - \overline x \right)^2 = \mathrm{SS}_{\mathrm W} + \mathrm{SS}_{\mathrm B}. \end{align}


D. Malchiodi, Superhero data science. Vol 1: probabilità e statistica: Analisi della varianza, 2017.
Powered by Jupyter Notebook