1 Variabili continue: grafici

Ricorda di scegliere il grafico in base al tipo di variabile:

  • barplot (grafico a barre) e pie plot (grafico a torta) per qualitative
  • istogramma per quantitative continue o discrete con un range elevato di valori
  • boxplot (diagramma a scatole e baffi) per il confronto tra gruppi.
Grafico Variabile Funzione principale
Diagramma a barre Qualitativa barplot()
Grafico a torta Qualtitativa pie()
Diagramma a punti Quantitativa stripchart()
Istogramma Quantitativa continua hist()
Boxplot Quantitativa (confronto gruppi) boxplot()

💡 Se si tratta di variabili discrete (spesso di conteggio) che assumono solo valori bassi, quindi con pochi valori distinti (esempio sono la variabile numero di fratelli, o il numero di sinistri con l’auto che l’assicurato denuncia in un anno) l’uso dei semplici grafici visti per le variabili categoriali si rivela spesso del tutto adeguato.

Consideriamo il data frame Cars93 del pacchetto MASS (nel seguito ci restringeremo alle analisi di Price e Length).

2 Il diagramma a punti: stripchart

Rappresenta su un grafico i singoli valori osservati (quindi non si perde alcun dettaglio). La funzione stripchart() del package graphics può essere usata a tale scopo.

library(MASS)
data("Cars93")

par(mfrow=c(1,2))
stripchart(Cars93$Price, pch=19, main ="Price", xlab = "Price",
           method="stack", cex=1.2)
stripchart(Cars93$Length, pch=18, col = "red", main ="Length", xlab = "Length",
           method="stack", cex=1.2)

par(mfrow=c(1,1))

Il parametro method=stack consente di sovrapporre dati osservati che risultano spesso avere il medesimo valore. Il default è invece quello di mostrare il valore osservato senza dare conto della eventuale molteplicità.

3 Istogramma: hist()

L’istogramma divide il range di una variabile continua in intervalli (classi) e rappresenta la frequenza di ciascuna classe con un rettangolo proporzionale in altezza (o area, in caso di densità).

Il suo processo di costruzione prevede:

  • la suddivisione del range dei dati in intervalli (classi o bins);
  • il conteggio delle osservazioni che ricadono in ciascun intervallo;
  • la rappresentazione di tali frequenze tramite rettangoli contigui.
  • L’altezza di ogni rettangolo è quindi determinata in modo tale da consentire una corrispondenza fra l’area del rettangolo (\(A_i\)) che insiste su $(a_i, a_{i+1})] e la frequenza relativa \(f_i\): ovvero \(fi \propto A_i\). Di solito si calcola l’altezza così che l’area di ogni rettangolo sia pari esattamente alla frequenza relativa, così che \(A_i = f_i\), e in tal caso l’area complessiva all’interno dei rettangoli che compongono il grafico risulterà pari ad 1. Le altezze dei rettangoli sono quindi poste pari alla frequenza relativa nell’intervallo diviso per l’ampiezza dell’intervallo stesso. Esse non rappresentano quindi le frequenze relative per gli intervalli in questione bensì la densità di frequenze relative (salvo nel caso in cui le basi dei i rettangoli abbiano tutti ampiezza pari a 1).
Intervallo Frequenza assoluta Frequenza relativa
\((a_1, a_2]\) \(n_1\) \(f_1\)
\((a_2, a_3]\) \(n_2\) \(f_2\)
\(\vdots\) \(\vdots\) \(\vdots\)
\((a_i, a_{i+1}]\) \(n_i\) \(f_i\)
\(\vdots\) \(\vdots\) \(\vdots\)
\((a_I, a_{I+1}]\) \(n_I\) \(f_I\)

3.1 Istogramma di Base

Gli istogrammi si ottengono con hist (nella versione di default, date \(n\) osservazioni viene stabilito il numero di classi \(I\) usando la regola di Sturges e quindi si calcola l’ampiezza di ogni classe dividenzo il range per il numero di classi; di default vengono rappresentate le frequenze assolute, poichè freq=T)

par(mfrow=c(1,2))
hist(Cars93$Price,
     main   = "Distribuzione del Prezzo",
     xlab   = "Prezzo (migliaia di $)",
     ylab   = "Frequenza assoluta",
     col    = "lightblue",
     border = "white")

hist(Cars93$Length,
     main   = "Distribuzione della Lunghezza",
     xlab   = "Lunghezza (in pollici)",
     ylab   = "Frequenza assoluta",
     col    = "lightblue",
     border = "white")

par(mfrow=c(1,1))

L’istogramma ci consente di:

  • individuare la forma della distribuzione (simmetria o asimmetria)
  • osservare la dispersione dei valori
  • rilevare eventuali valori anomali (outlier)

💡 La distribuzione di Price è fortemente asimmetrica a destra: pochi valori molto elevati allungano la coda. Questo può essere tipico di variabili economiche.

Esempi di Asimmetria/Simmetria

3.2 Numero di Classi: il Parametro breaks

Il numero di classi influenza notevolmente l’aspetto dell’istogramma. Troppo poche classi nascondono la forma della distribuzione; troppo la frammentano:

library(insuranceData)
data("AutoBi")
par(mfrow = c(1, 2))

hist(Cars93$Price,
     breaks = 5,
     main   = "breaks = 5",
     xlab   = "Prezzo (migliaia di $)",
     ylab   = "Frequenza assoluta",
     col    = "lightblue",
     border = "white")

hist(Cars93$Price,
     breaks = 50,
     main   = "breaks = 50",
     xlab   = "Prezzo (migliaia di $)",
     ylab   = "Frequenza assoluta",
     col    = "lightblue",
     border = "white")

par(mfrow = c(1, 1))

3.3 Densità e Curva Sovrapposta

Con freq = FALSE sull’asse delle ordinate compare la densità (area totale = 1). Si può sovrapporre una curva di densità stimata con density():

  • fornisce una rappresentazione continua e smussata della distribuzione;
  • non dipende dalla suddivisione in classi (meno sensibile ai breaks);
  • aiuta a cogliere meglio la forma generale dei dati.
hist(Cars93$Price,
     freq = FALSE, 
     breaks = 20,
     main   = "Istogramma e densità",
     xlab   = "Prezzo (migliaia di $)",
     ylab   = "Frequenza assoluta",
     col    = "lightblue",
     border = "white")
lines(density(Cars93$Price, na.rm = TRUE),
      col = "darkblue",
      lwd = 2)

3.4 Trasformazione Logaritmica

Quando la distribuzione è molto asimmetrica, la trasformazione logaritmica rende la forma più leggibile:

hist(log(Cars93$Price),
     freq   = FALSE,
     main   = "Distribuzione di log(Price)",
     xlab   = "log(Price)",
     col    = "lightblue",
     border = "white")
lines(density(log(Cars93$Price)),
      col = "darkblue",
      lwd = 2)

3.5 Classi di Ampiezza Diversa

Se le classi hanno ampiezza diversa, è obbligatorio usare freq = FALSE (densità): con frequenze assolute, le classi più ampie apparirebbero artificialmente più alte.

breaks_custom <- c(0, 10, 30,  max(Cars93$Price))

hist(Cars93$Price,
     breaks = breaks_custom,
     freq   = FALSE,
     main   = "Istogramma con classi di ampiezza diversa",
     xlab   = "Price",
     col    = "lightblue",
     border = "white")

⚠️ Con classi di ampiezza uguale freq = TRUE e freq = FALSE danno istogrammi con la stessa forma (cambia solo la scala dell’asse y). Con classi di ampiezza diversa, usare sempre freq = FALSE per non distorcere la rappresentazione.

3.5.1 Esercizio

Esercizio 1

Usando la variabile LOSS di AutoBi:

  1. Traccia l’istogramma con il numero di classi di default.

  2. Traccia l’istogramma con breaks = 20 e con breaks = 5 affiancati (par(mfrow = ...)).

  3. Traccia l’istogramma di log(LOSS) con la curva di densità sovrapposta.

  4. La distribuzione di LOSS è simmetrica? Come cambia dopo la trasformazione log?

4 Boxplot: boxplot()

Il boxplot (diagramma a scatola e baffi) riassume graficamente cinque statistiche:

         ┌─────────────┐
  ──┬────┤     IQR     ├────┬──  ×  ×  (outlier)
    │    └─────────────┘    │
  Q1-1.5·IQR  Q1  Me  Q3  Q3+1.5·IQR
Elemento Statistica
Linea centrale Mediana (Q2)
Bordi della scatola Primo (Q1) e terzo quartile (Q3)
Baffi Fino a 1.5 × IQR dal bordo della scatola
Punti isolati Outlier (oltre i baffi)

La definizione di outlier nel diagramma a scatola deriva dalla seguente regola:

  • sono outlier (valori anomali) quei valori che sono più distanti dai bordi della scatola (cioè dai quartili) più di una volta e mezza lo scarto interquartile IQR. Pertanto tutti i punti che sono superiori a Q3 + 1.5IQR o inferiori a Q1 - 1.5IQR verranno annotati separatamente sul grafico;

  • se essi esistono (a destra o a sinistra) di conseguenza il baffo non va esteso fino al massimo o al minimo valore osservato, va invece esteso:

    – a destra, fino al più grande valore che non sia segnalato come outlier (cioè il massimo valore osservato che risulti inferiore a Q3 + 1.5IQR);

    – a sinistra, fino al più piccolo valore che non sia segnalato come outlier (cioè il minimo valore osservato che risulti superiore a Q1 - 1.5IQR).

4.1 Quantili e Medie di Posizione

Le medie di posizione sono valori che hanno una posizione definita rispetto alla sequenza ordinata dei dati. Il concetto chiave è quello di quantile empirico.

Fissata una proporzione \(p\) (con \(0 \leq p \leq 1\)), si definisce \(x_p\) quantile empirico di ordine \(p\) come quel valore \(x \in \mathbb{R}\) tale che:

\[x_p : \frac{\#(x_i \leq x_p)}{n} = p\]

⚠️ Non è detto che esista un unico valore che realizzi la condizione, né che ne esista alcuno. In R si usa l’interpolazione lineare fra i due quantili successivi osservati. Si veda ?quantile per le diverse convenzioni disponibili.

4.1.1 La Mediana e i Quartili

I quantili \(x_{0.25}\), \(x_{0.5}\), \(x_{0.75}\) sono detti rispettivamente primo, secondo e terzo quartile. Il secondo quartile è anche detto mediana (\(Me\)).

La mediana è definita convenzionalmente come:

\[Me = \begin{cases} x_{\left(\frac{n+1}{2}\right)} & \text{se } n \text{ è dispari} \\ \dfrac{x_{\left(\frac{n}{2}\right)} + x_{\left(\frac{n}{2}+1\right)}}{2} & \text{se } n \text{ è pari} \end{cases}\]

La mediana (a differenza) dell a media è molto meno sensibile ai valori estremi.

⚠️ Si ricorda che come per tutte le funzioni statistiche, in presenza di NA occorre usare na.rm = TRUE

median(Cars93$Price,   na.rm = TRUE)
#> [1] 17.7
median(Cars93$Length, na.rm = TRUE)
#> [1] 183
# Quartili (e quantili di ordine 0 e 1, cioè minimo e massimo)
quantile(Cars93$Price, na.rm = TRUE)
#>   0%  25%  50%  75% 100% 
#>  7.4 12.2 17.7 23.3 61.9
# Decili
quantile(Cars93$Price, probs = seq(0, 1, 0.1), 
         na.rm = TRUE, digits = 2)
#>    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
#>  7.40  9.84 11.34 13.74 15.78 17.70 19.34 20.84 26.22 33.62 61.90
# Verifica: il quantile di ordine 0.5 è la mediana
median(Cars93$Price) == quantile(Cars93$Price, probs = 0.5)
#>  50% 
#> TRUE

4.1.2 Scarto Interquartile (IQR)

Lo scarto interquartile è una misura di variabilità data dalla differenza fra il terzo e il primo quartile: \[IQR = Q_3 - Q_1\]

Tiene conto solo della dispersione nella parte centrale della distribuzione (il 50% dei dati centrali), ignorando i valori estremi. È leggibile direttamente dall’output di summary() ma esiste una funzione dedicata IQR() per il suo calcolo.

IQR(Cars93$Price,   na.rm = TRUE)
#> [1] 11.1
IQR(Cars93$Length, na.rm = TRUE)
#> [1] 18

4.2 Boxplot singolo in R

La funzione per ottenere un boxplot in R è boxplot().

boxplot(Cars93$Price,
        main    = "Boxplot del  Prezzo",
        ylab    = "Prezzo (migliaia di $)",
        col     = "lightblue")

4.3 Boxplot Multipli: Confronto tra Gruppi

L’uso più efficace del boxplot è il confronto della distribuzione di una variabile quantitativa tra i gruppi definiti da una variabile qualitativa. Si usa la notazione formula y ~ gruppo e l’argomento data permette di evitare l’accesso mediante dollaro alle variabili. L’argomento las = 2 chepermette di ruotare le etichette sull’asse x (utile quando sono lunghe)

boxplot(Price ~ Type,
        data = Cars93,
        main = "Prezzi per tipo di auto",
        xlab = "Tipo",
        ylab = "Prezzi (migliaia di $)", 
        las = 2)

L’argomento notch = TRUE aggiunge le incisioni per il confronto visivo delle mediane.

boxplot(Price ~ Origin,
        data  = Cars93,
        main  = "Prezzi per origine",
        xlab  = "Origine",
        ylab  = "Prezzi",
        col   = c("lightblue", "lightpink"),
        notch = TRUE)

💡 Se le incisioni di due boxplot non si sovrappongono, c’è evidenza visiva che le mediane siano diverse. In questo caso non vi è evidenza visiva che le mediane siano diverse.

4.3.1 Esercizio

Esercizio 2

Usando AutoBi:

  1. Costruisci i boxplot di LOSS e log(LOSS)

  2. Costruisci i boxplot di log(LOSS) separatamente per Ricorso all’avvocato (ATTORNEY). Ci sono differenze tra i danni nei due gruppi? Questo conferma quanto visto rendendo LOSS in factor

  3. Costruisci i boxplot di CLMAGE per le combinazioni di sesso e avvocato usando interaction(CLMSEX, ATTORNEY).

5 Funzione di ripartizione empirica: ecdf()

Una modo per rappresentare graficamente un insieme di dati è quello di ottenere una versione empirica della funzione di ripartizione di una variabile aleatroria \(X\), cioè \(F(x) = Pr(X \leq x)\). Inoltre l’inversa della funzione di ripartizione fornisce i quantili per cui il quantile \(x_p = F^{-1}(p)\). Ovviamente il quantile per ogni \(0 < p < 1\) esiste se la funzione di ripartizione non ha discontinuità.

L’equivalente empirico di tale funzione, denotato con \(F_n(x)\), avendo osservato l’insieme di dati \(x_1, x_2, \ldots , x_n\), per una variabile quantitativa, si ottiene cercando il valore che \(x \in \mathbb{R}\) fornisce la proporzione di unità inferiori o pari a \(x\), ovvero

\[F_n(x) = \text{proporzione}(x_i \leq x) = \frac{\text{numero di valori} \leq x}{n}\]

Considerando i dati ordinati \(x_{(1)} < x_{(2)} < \ldots < x_{(n-1)} < x_{(n)}\) la funzione di ripartizione empirica (cumulata empirica) sarà pari a:

  • \(F_n(x) = 0\) se \(x \leq x_{(1)}\);
  • in \(x = x_{(1)}\) la funzione fa un salto ed è pari a \(1/n\) e rimane costante fino al prossimo valore \(x_{(2)}\);
  • in \(x = x_{(2)}\) la funzione fa un ulteriore salto ancora di altezza \(1/n\);
  • e così via fino a \(x_{(n)}\) partire dal quale la funzione assume il valore 1;
  • se vi sono valori ripetuti il salto sarà pari a \(m/n\) ove \(m\) è il numero di valori ripetuti.

Si tratta quindi di una funzione a scalini definita come: \[F_n(x) = \begin{cases} 0 \quad \text{se} \quad x < x(1) \\ i/n \quad \text{se} \quad x_{(i)} \leq x < x(i+1) \quad i = 1, 2, \ldots , n - 1 \\ 1 \quad \text{se} \quad x \geq x_{(n)} \end{cases}\]

In R otteniamo la funzione di ripartizione empirica mediante ecdf() che è basata sulla funzione stepfun() che è appositamente costruita per trattare funzioni a gradini.

Esempio

par(mfrow=c(1,2))
# consideriamo l!insieme di 9 dati nel vettore xx
xx<-c(148, 172, 155,  168, 162,  175, 178, 186, 174)
sort(xx)
#> [1] 148 155 162 168 172 174 175 178 186
plot(ecdf(xx))

# Aggiungiamo un valore ripetuto
yy <- c(xx, 168)
plot(ecdf(yy))

Nel caso della lunghezza delle auto, si ha

plot(ecdf(Cars93$Length), main="funzione di ripartizione empirica per lunghezza",
xlab="lunghezza", ylab="proporzione di casi")

Tale rappresentazione grafica è di interesse per vari motivi:

  • conserva l’informazione su ogni singolo valore;
  • è del tutto adeguata anche con variabili quantitative discrete;
  • permette di visualizzare dove si trovi la mediana o alcuni percentili (la funzione dei quantili è in effetti l’inversa della funzione di ripartizione). Da questo si intende chiaramente come solo alcuni quantili sono definiti, mentre per altri occorre ricorrere ad approssimazioni (ad esempio considerando interpolazioni lineari fra i due quantili con l’ordine più prossimo);
  • la sua interpretazione non è agevole e immediata, ma è evidente che i punti in cui la curva cresce più veleocemente sono quelli in cui si sono osservati molti casi;
  • è estremamente vantaggiosa se si vuole confrontare la distribuzione empirica con una distribuzione teorica (aspetto che verrà trattato più avanti).
plot(ecdf(Cars93$Length), main="funzione di ripartizione empirica per lunghezza e mediana",
xlab="lunghezza", ylab="proporzione di casi")
segments(130, 0.5,  
         median(Cars93$Length), 0.5,  
        col = "red", lty = "dashed")
segments(median(Cars93$Length), 0,
         median(Cars93$Length), 0.5,
        col = "red", lty = "dashed")