Funzione di distribuzione empirica discreta

La funzione viene definita a partire dalle frequenze relative cumulate.
Prendiamo in considerazione una delle nostre variabili quantitative come ad esempio la percentuale di individui che hanno conseguito una laurea per ogni regione. Diciamo che i valori distinti che può assumere questa variabile sono z1,z2 … zk e assumiamo che siano ordinati in ordine crescente. Consideriamo poi le frequenze relative e le frequenze relative cumulate:

dove la generica \(F_i\) rappresenta la proporzione dei dati del campione minori o uguali di zi. Una funzione di distribuzione empirica discreta risulta:

La funzione è definita per ogni x reale e ha diverse caratteristiche:

• È una funzione a scalini non decrescente
• In corrispondenza di ogni salto assume il valore a sinistra sull’asse delle ascisse
• La funzione vale:
  • 0 per ogni valore minore dell’osservazione minima
  • 1 per ogni valore maggiore dell’osservazione massima

R mette a disposizione la funzione ecdf() (empirical cumulative distribution function) che ci permette di disegnare il grafico della funzione e di determinare per ogni x reale il valore di tale funzione.

Dato che i valori contenuti nei nostri vettori sono tutti diversi, la funzione non ci è molto di aiuto nella nostra indagine.
Vediamo un esempio di applicazione sul vettore relativo ai laureati per regione in Italia:

plot(
  ecdf(laurea),
  main="Funzione di distribuzione empirica discreta relativo alle lauree ",
  verticals = TRUE , col="red")

Nel caso si voglia sapere il valore della funzione in un dato punto x, R permette di farlo nel seguente modo:

ecdf(laurea)(17)

## [1] 0.3

Funzione di distribuzione empirica continua

La funzione di distribuzione empirica continua viene invece utilizzata quando si lavora con dati raccolti in classi. Infatti considerando la tipologia del nostro campione di dati è più opportuno suddividere le informazioni a disposizione in k distinte classi. Anche in questo caso la funzione viene definita a partire dalle frequenze relative cumulate.

La funzione di distribuzione empirica continua è dunque così definita: La funzione è definita per ogni x reale e:

• Vale 0 per ogni x minore di z1
• Vale 1 per ogni x maggiore o uguale a zk+1
• Se zi < x < zi+1 coincide con il segmento che passa per i punti (zi, Fi) e (zi+1,Fi+1)

Introduciamo dunque le seguenti classi:
[0,5) [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) [35,40) [40,45) [45,50)

classiDistrEmpCont <-c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50)

Descriviamo ora come ottenere il grafico della funzione di distribuzione empirica continua per quanto riguarda la licenza elementare:

FrequenzaCumulateLicEle <- cumsum (table (
  cut(licenza_elementare_o_nessuno ,
  breaks =classiDistrEmpCont ,
  right =FALSE )))/length (licenza_elementare_o_nessuno)

FrequenzaCumulateLicEle <- c(0,FrequenzaCumulateLicEle)
plot(classiDistrEmpCont , 
     FrequenzaCumulateLicEle , type = "b", axes = FALSE , 
     main = "Funzione di distribuzione empirica continua licenza elementare",
     col="red", xlab = "Classi", ylab = "Frequenze cumulate")
axis(1, classiDistrEmpCont, cex.axis=0.80)
axis(2, format (FrequenzaCumulateLicEle , digits = 2), cex.axis=0.80)
box ()

È possibile notare dal grafico come ci sia una forte concentrazione dei dati, il 65% dei dati è nelle prime due classi, e l’unità si raggiunge velocemente in quelle successive.

Procediamo con i grafici per ogni altra variabile.

I dati sono concentrati principalmente verso le ultime classi, i primi valori si riscontrano nella classe [25,30). In generale abbiamo percentuali mediamente molto alte.

Analogamente alla prima situazione anche qui l’unità si raggiunge subito e c’è una concentrazione di dati nelle prime classi: addirittura con le prime due classi copriamo l’80% dei dati.

I dati sono concentrati principalmente verso le ultime classi, le percentuali sono mediamente molto alte.

Il grafico si differenzia dagli altri per raggiungere l’unità nelle classi centrali, i dati sono molto più centrati rispetto agli altri dove abbiamo visto una concentrazione nelle classi iniziali e quelle finali, qui invece la concentrazione è in quelle centrali.

Media campionaria

Supponiamo di avere un campione di numerosità N con valori x1,x2 … xn.
Definiamo la media campionaria come la media aritmetica di questi valori:

La media campionaria ha diverse proprietà:

• Gode della proprietà di linearità:

• È una media pesata dei valori distinti assunti dai dati dove ogni valore distinto usa come peso la frequenza:

• La somma algebrica degli scarti della media campionaria è sempre nulla:

• La media è influenzata da tutti i dati e in particolar modo da valori particolarmente grandi o piccoli, cioè da valori anomali

R ci mette a disposizione la funzione mean() per calcolare la media campionaria dei nostri dati.
Applichiamo questa funzione ad ogni vettore che compone il nostro campione:

mediaCampionariaElementari <- mean(licenza_elementare_o_nessuno)
mediaCampionariaElementari

## [1] 5.925

mediaCampionariaMedie <- mean(licenza_media)
mediaCampionariaMedie

## [1] 33.115

mediaCampionariaDiploma <- mean(diploma)
mediaCampionariaDiploma

## [1] 6.755

mediaCampionariaMaturita <- mean(maturita)
mediaCampionariaMaturita

## [1] 35.68

mediaCampionariaLaurea <- mean(laurea)
mediaCampionariaLaurea

## [1] 18.51

Mediana campionaria

Un altro indice di posizione è la mediana.
Data una distribuzione di un carattere quantitativo oppure qualitativo ordinabile si definisce mediana (o valore mediano) il valore che separa la metà più grande e la metà più piccola di un campione di dati, cioè il valore che si trova nel mezzo della distribuzione.

I passi per calcolare la mediana sono i seguenti:

• Si ordinano gli n elementi in ordine crescente
• Se il numero di dati è dispari la mediana corrisponde al valore \(\frac{n+1}{2}\)
• Se il numero di dati è pari allora la mediana è calcolato come la media aritmetica dei valori alla posizione \(\frac{n}{2}\) e alla posizione \(\frac{n}{2}+1\)

È importante notare come la mediana campionaria dipenda solo da uno o due valori e quindi non risente degli estremi.

In R per calcolare la mediana abbiamo a disposizione la funzione median().
Calcoliamola dunque per i nostri vettori:

medianaElementari <- median(licenza_elementare_o_nessuno)
medianaElementari

## [1] 4.45

medianaMedie <- median(licenza_media)
medianaMedie

## [1] 31.95

medianaDiploma <- median(diploma)
medianaDiploma

## [1] 5.4

medianaMaturita <- median(maturita)
medianaMaturita

## [1] 36.05

medianaLaurea <- median(laurea)
medianaLaurea

## [1] 18.55

Visualizziamo in forma tabellare i risultati:

Dalla tabella notiamo che il calcolo della media e della mediana ha prodotto risultati simili per ciascuna categoria, ad esempio il vettore laurea dà praticamente lo stesso risultato per media e mediana, mentre per la prima variabile vediamo la differenza maggiore.
Confrontare media e mediana è utile poiché se queste misure sono simili, come nel nostro caso, la distribuzione di frequenza tende ad essere simmetrica.
Vedremo nel paragrafo in cui studiamo la forma della distribuzione se questa intuizione sarà verificata.

Moda campionaria

La moda identifica la modalità che si presenta con la maggiore frequenza (assoluta o relativa) in un campione. Se ci sono più modalità con frequenza massima, ciascuna viene definita come valore modale.
Sulla moda c’è da dire che questa, a differenza di media e mediana, può essere calcolata anche quando si ha a che fare con dati di tipo qualitativo.
Nel nostro caso è inutile calcolarla in quanto abbiamo a che fare con percentuali.

Infine sottolineiamo che media, mediana e moda sono detti indici di posizione centrale poiché descrivono attorno a quali valori è centrato l’insieme di dati.

Quantili

Sono definiti come indici di sintesi di posizione non centrali.
I quantili dividono l’insieme dei dati in un fissato numero di parti uguali.
Definiamo i quantili in base a come viene suddiviso l’insieme dei dati:

• Se l’insieme dei dati viene suddiviso in 100 parti mediante 99 quantili si parla di percentili
• Se l’insieme dei dati viene suddiviso in 10 parti mediante 9 quantili si parla di decili
• Se l’insieme dei dati viene suddiviso in 4 parti mediante 3 quantili si parla di quartili

Ovviamente sia quartili che decili sono un caso particolare dei percentili.

In R esistono 9 algoritmi per calcolare i quantili utilizzando la funzione quantile() e specificando il parametro type con un numero nell’insieme {1 … 9}. R utilizza di default l’algoritmo di tipo 7. Ad ogni modo se il numero di osservazioni è elevato, i valori calcolati dagli algoritmi tendono a coincidere.

I percentili maggiormente utilizzati sono il 25-esimo, il 50-esimo e il 75-esimo detti rispettivamente primo quartile (Q1), secondo quartile (Q2), terzo quartile (Q3).

Vediamo dunque il calcolo dei quartili per il vettore laurea come cambia anche in base al tipo di algoritmo scelto:

quantile(laurea)

##     0%    25%    50%    75%   100% 
## 13.600 16.425 18.550 20.475 24.900

quantile(laurea, type = 2)

##    0%   25%   50%   75%  100% 
## 13.60 16.35 18.55 20.55 24.90

quantile(laurea, type = 1)

##   0%  25%  50%  75% 100% 
## 13.6 16.2 18.3 20.4 24.9

Varianza e deviazione standard campionaria

Gli indici di posizione sono importanti, ma le informazioni che ci danno non tengono conto della variabilità dei dati: può capitare di avere a che fare con distribuzioni di frequenza molto diverse tra di loro ma che presentano la stessa media campionaria ad esempio.
È dunque importante arricchire l’indagine su un dato campione di dati con degli indici per capaci di misurare la variabilità dei dati.

Introduciamo dunque la varianza campionaria e la deviazione standard campionaria.

La varianza campionaria è indicata con \(s^2\) ed è definita:

La deviazione standard è invece indicata con \(s\) ed è definita come la radice quadrata della varianza campionaria:

Varianza e deviazione standard sono detti indici di dispersione in quanto misurano la dispersione dei dati intorno alla media.
La deviazione standard ha la stessa unità di misura dei valori osservati, mentre la varianza ha come unità di misura il quadrato dell’unità di misura dei valori di riferimento.

La deviazione standard viene detta anche scarto quadratico medio per come è definita.

Le proprietà della varianza sono le seguenti:

• Non gode della proprietà di linearità
• Sommare una costante a ciascuno dei dati non fa cambiare la varianza
• Moltiplicare i dati del campione per un fattore costante fa sì che la varianza campionaria risulti moltiplicata per il quadrato di tale fattore

Il valore della varianza è 0 quando tutti i valori sono uguali tra loro e risulta essere più grande tanto più le osservazioni si discostano dalla media.

In R per calcolare la varianza campionaria usiamo il comando var(), mentre per calcolare la deviazione standard usiamo la funzione sd().
Valutiamo quindi quanto si disperdono i nostri dati rispetto alla media.
Calcoliamo prima la varianza per ogni categoria:

varianzaElementari <- var(licenza_elementare_o_nessuno)
varianzaElementari

## [1] 8.2525

varianzaMedie <- var(licenza_media)
varianzaMedie

## [1] 18.04134

varianzaDiploma <- var(diploma)
varianzaDiploma

## [1] 19.25734

varianzaMaturita <- var(maturita)
varianzaMaturita

## [1] 11.42695

varianzaLaurea <- var(laurea)
varianzaLaurea

## [1] 7.798842

Calcoliamo adesso la deviazione standard:

deviazioneStandardElementari <- sd(licenza_elementare_o_nessuno)
deviazioneStandardElementari

## [1] 2.872716

deviazioneStandardMedie <- sd(licenza_media)
deviazioneStandardMedie

## [1] 4.24751

deviazioneStandardDiploma <- sd(diploma)
deviazioneStandardDiploma

## [1] 4.388319

deviazioneStandardMaturita <- sd(maturita)
deviazioneStandardMaturita

## [1] 3.380377

deviazioneStandardLaurea <- sd(laurea)
deviazioneStandardLaurea

## [1] 2.792641

Questi indici di dispersione sono relativi ai singoli dati presi in considerazione, se vogliamo però confrontare le variazioni tra diversi insiemi di dati si utilizza il coefficiente di variazione.

Coefficiente di variazione

Il coefficiente di variazione è definito come il rapporto tra la deviazione standard campionaria e il modulo della media campionaria:

Il coefficiente di variazione è un numero puro, è un indice adimensionale, cioè non dipende dall’unità di misura utilizzata (la media campionaria e la deviazione standard campionaria sono espressi in identiche unità di misura).
È facile capire che il coefficiente di variazione ha senso che sia calcolato solo con campioni con media campionaria non nulla.
Non dobbiamo guardare alla varianza in senso assoluto, ma dobbiamo confrontarla con la media per fare una buona valutazione.

Nel nostro caso è importante guardare al coefficiente di variazione perché ci interessa confrontare insiemi che hanno differenti range di variazione, cioè insiemi in cui la differenza tra il massimo e il minimo è molto diversa.

Il coefficiente risulta essere sempre maggiore di 0, ma se compreso tra 0 e 1 vuol dire che la media è più grande della deviazione standard.
Invece più è grande più c’è dispersione nei dati e quindi il valore medio non è molto significativo.

Dato che in R non c’è una funzione per calcolare il coefficiente di variazione, introduciamo la seguente funzione:

cv <- function (x){
        sd(x)/abs(mean(x))
      }

Utilizziamo quindi la funzione per ricavare l’indice per i nostri vettori:

cvElementare <- cv(licenza_elementare_o_nessuno)
cvMedie <- cv(licenza_media)
cvDiploma <- cv(diploma)
cvMaturita <- cv(maturita)
cvLaurea <- cv(laurea)

Utilizziamo quindi la seguente tabella per avere un quadro totale degli indici ricavati per fare le nostre opportune analisi, in ordine di riga sono riportate media, varianza, deviazione standard e coefficiente di variazione:

Dalla tabella notiamo come il coefficiente di variazione è minore per le variabili relative ai titoli medie, maturità e laurea, mentre è decisamente più alto per nessun titolo e soprattutto diploma dove si nota la variazione più grande tra i dati.

Skewness

Dato un insieme di dati numerici, si definisce skewness:

dove \(m_j\) è il momento centrato campionario definito:

In base al valore della skewness abbiamo che:

• Se \(\gamma1 = 0\), la distribuzione di frequenza è simmetrica
• Se \(\gamma1 > 0\), la distribuzione di frequenza ha la coda di destra più allungata per l’asimmetria positiva
• Se \(\gamma1 < 0\), la distribuzione di frequenza ha la coda di sinistra più allungata per l’asimmetria negativa

Definiamo una funzione per calcolare l’indice:

skw <-function (x){
  n<-length (x)
  m2 <-(n -1)*var(x)/n
  m3 <- (sum( (x-mean(x))^3))/n
  m3/(m2 ^1.5)
}

E calcoliamo i valori:

skw(licenza_elementare_o_nessuno)

## [1] 1.062309

skw(licenza_media)

## [1] 0.8847151

skw(diploma)

## [1] 1.268479

skw(maturita)

## [1] 0.020537

skw(laurea)

## [1] 0.08639237

Notiamo che l’indice non è mai pari a 0, quindi non c’è mai simmetria.
Se guardiamo attentamente i risultati notiamo come il vettore maturità mantenga la simmetria maggiore avvicinandosi allo 0, mentre tutti gli altri hanno una distribuzione che presenta una forte asimmetria positiva e quindi la coda è allungata verso destra.

Curtosi

L’indice di curtosi invece permette di misurare la densità dei dati intorno alla media ed è definita:

dove \(\beta2\) è :

Quello che facciamo tramite l’indice di curtosi è confrontare la distribuzione di frequenza dei nostri dati con una densità di probabilità normale caratterizzata da \(\beta_2 = 3\) e indice di curtosi \(\gamma_2 = 0\).

Risulta quindi:

• \(\beta_2 < 3\), \(\gamma_2 < 0\): la distribuzione di frequenza si definisce platicurtica, ossia la distribuzione di frequenza è più piatta di una normale
• \(\beta_2 > 3\), \(\gamma_2 > 0\): la distribuzione di frequenza si definisce leptocurtica, ossia la distribuzione di frequenza è più piccata di una normale
• \(\beta_2 = 3\), \(\gamma_2 = 0\): la distribuzione di frequenza si definisce normocurtica, ossia segue la curva di una normale

Definiamo quindi la funzione per calcolare l’indice nel seguente modo:

curt <-function (x){
  n<-length (x)
  m2 <-(n-1)*var (x)/n
  m4 <- (sum( (x-mean(x))^4) )/n
  m4/(m2^2) -3
}

E calcoliamo i valori:

curt(licenza_elementare_o_nessuno)

## [1] -0.4591248

curt(licenza_media)

## [1] 0.09400184

curt(diploma)

## [1] 1.492363

curt(maturita)

## [1] -0.7978636

curt(laurea)

## [1] -0.1815536

La curtosi risulta quindi negativa per la prima, la quarta e la quinta categoria, cioè platicurtica. Negli altri due casi è invece leptocurtica, anche se c’è da notare che il vettore licenza media si avvicina ad avere una forma normale.

Residui

Dopo aver trovato la retta di regressione, possiamo osservare qual è il discostamento tra i valori osservati (le coppie \((x_i,y_i)\)) e i valori stimati (le coppie \((x_i,\hat{y_i})\)).

I valori stimati sono espressi secondo l’equazione \[\hat{y_i} = \alpha + \beta x_i\] ottenuti mediante la retta di regressione. Risulta inoltre che la media campionaria dei valori stimati è uguale alla media campionaria dei valori osservati.

I residui sono dunque definiti come \[E_i = y_i - \hat{y_i} = y_i - (\alpha + \beta x_i)\] e mostrano di quanto si discostano valori osservati e valori stimati.

Si nota sui residui:

• La media campionaria dei residui risulta sempre nulla
• La varianza dei residui è \(S_E^2 = \frac{1}{n-1}\sum_{i=1}^n E_i^2\) in quanto la media campionaria è 0

Per calcolare il vettore dei valori stimati utilizziamo la funzione fitted(),passando come argomento lm(y~x):

valoriStimati <- fitted(lm(laurea~licenza_media))
valoriStimati

##        1        2        3        4        5        6        7        8 
## 17.63255 16.95558 19.87174 19.19478 20.54871 18.98648 20.65286 20.28834 
##        9       10       11       12       13       14       15       16 
## 19.61137 21.01738 19.35100 21.53812 19.50722 18.15329 16.79936 15.18506 
##       17       18       19       20 
## 19.03855 17.84085 14.92469 13.10209

I valori stimati possono essere ottenuti anche attraverso l’attributo “fitted.values” di linear model:

linearModel$fitted.values

##        1        2        3        4        5        6        7        8 
## 17.63255 16.95558 19.87174 19.19478 20.54871 18.98648 20.65286 20.28834 
##        9       10       11       12       13       14       15       16 
## 19.61137 21.01738 19.35100 21.53812 19.50722 18.15329 16.79936 15.18506 
##       17       18       19       20 
## 19.03855 17.84085 14.92469 13.10209

Per calcolare i residui invece di usa la funzione resid() passando anche in questo caso lo stesso argomento:

valoriResidui <- resid(lm(laurea~licenza_media))
valoriResidui

##          1          2          3          4          5          6 
##  0.4674512  1.0444166 -0.2717418  1.2052235 -1.7487071 -1.0864794 
##          7          8          9         10         11         12 
## -0.8528556  0.9116641  1.0886295  0.1826246  1.3490008  3.3618820 
##         13         14         15         16         17         18 
## -1.2072220  1.5467087 -1.2993607 -1.1850587 -2.5385537 -2.7408458 
##         19         20 
## -1.3246875  3.0979115

I residui possono essere ottenuti anche attraverso l’attributo “residuals” di linear model:

linearModel$residuals

##          1          2          3          4          5          6 
##  0.4674512  1.0444166 -0.2717418  1.2052235 -1.7487071 -1.0864794 
##          7          8          9         10         11         12 
## -0.8528556  0.9116641  1.0886295  0.1826246  1.3490008  3.3618820 
##         13         14         15         16         17         18 
## -1.2072220  1.5467087 -1.2993607 -1.1850587 -2.5385537 -2.7408458 
##         19         20 
## -1.3246875  3.0979115

Dei residui possiamo calcolare mediana, varianza e deviazione standard, mentre non è possibile calcolare il coefficiente di variazione in quanto la media dei residui è 0.
Vediamo dunque nel nostro caso:

median(linearModel$residuals)

## [1] -0.04455862

var(linearModel$residuals)

## [1] 2.906521

sd(linearModel$residuals)

## [1] 1.704852

Rappresentiamo dunque graficamente i valori dei residui, vediamo tre possibili rappresentazioni:

1. aggiungiamo al grafico dello scatterpolot e la retta di regressione dei segmenti verticali che visualizzano i residui: congiungiamo i punti \((x_i,y_i)\)) e i punti \((x_i,\hat{y_i})\)

2. Usiamo un diagramma dei residui: un grafico in cui i valori dei residui sono posti sull’asse delle ordinate e quelli della variabile indipendente sull’asse delle ascisse Il diagramma dei residui aiuta a comprendere quale è l’adattamento della retta di regressione rispetto ai dati, consentendo di identificare quali sono le informazioni che hanno una forte influenza sulla collocazione e direzione della retta di regressione.
   Occorre notare che la posizione della retta di regressione è fortemente influenzata dalla presenza di eventuali valori anomali che si discostano in modo significativo dagli altri. L’analisi dei residui aiuta ad individuare eventuali punti isolati (valori anomali) dovuti ad errori nella stima. Tali valori possono perturbare significativamente la stima dei parametri di regressione e influenzare l’interpretazione dei residui. Eliminando i valori anomali la varianza campionaria dei residui diminuisce.

3. Rappresentare i residui standardizzati rispetto ai valori stimati. I residui standardizzati sono definiti: \[E_i^{(s)} = \frac{E_i - \overline{E}}{s_E} = \frac{E_i}{s_E}\] caratterizzati da media nulla e deviazione standard pari a 1.

I punti indicano la posizione dove si collocano i residui standardizzati rispetto ai valori stimati con la retta di regressione. La retta orizzontale è posizionata nello zero, che corrisponde alla media campionaria dei residui standardizzati. Anche in questo caso i punti sono disposti quasi casualmente attorno alla linea orizzontale e non si evidenzia nessuna tendenza particolare nella distribuzione dei punti.

Coefficiente di determinazione

Un altro indice da considerare è il coefficiente di determinazione, definito come:

È definita cioè come il rapporto tra varianza dei valori stimati (con la retta di regressione) e varianza dei valori osservati.

Il coefficiente varia tra 0 e 1:

• Vicino ad 1 significa che i punti tenderanno ad allinearsi lungo la retta di regressione
• Vicino a 0 non c’è una retta capace di interpolare i punti

Nel caso della regressione lineare semplice il coefficiente di determinazione è il quadrato del coefficiente di correlazione, cioè \[D^2 = r_{xy}^2\]

Per calcolarlo in R possiamo procedere nei due modi seguenti:

(cor(licenza_media,laurea))^2

## [1] 0.6273138

summary(linearModel)$r.square

## [1] 0.6273138

Residui

I residui come detto mostrano di quanto si discostano i valori osservati da quelli stimati con la retta di regressione.
Sono definiti come \[E_i = y_i - \hat{y_i} = y_i - (\alpha + \beta_1 x_{i,1}+ \beta_2 x_{i,2}+...+ \beta_p x_{i,p})\]

In modo analogo con la regressione lineare semplice, si nota che:

• La media campionaria risulta sempre nulla
• La varianza dei residui è \(S_E^2 = \frac{1}{n-1}\sum_{i=1}^n E_i^2\) dato che la media campionaria è 0

Per calcolare il vettore dei valori stimati utilizziamo la funzione fitted(),passando come argomento lm(y~x1 + x2 … +xp):

valoriStimatiMult <- fitted(linearModelMultiplo)
valoriStimatiMult

##        1        2        3        4        5        6        7        8 
## 18.08645 17.98388 19.59623 20.40465 18.90495 17.78806 19.80304 21.11115 
##        9       10       11       12       13       14       15       16 
## 20.70246 21.20794 20.70477 24.93436 18.28509 19.68865 15.45004 13.93727 
##       17       18       19       20 
## 16.56616 15.14733 13.63613 16.26136

I valori stimati possono essere ottenuti anche attraverso l’attributo “fitted.values” di linear model:

linearModelMultiplo$fitted.values

##        1        2        3        4        5        6        7        8 
## 18.08645 17.98388 19.59623 20.40465 18.90495 17.78806 19.80304 21.11115 
##        9       10       11       12       13       14       15       16 
## 20.70246 21.20794 20.70477 24.93436 18.28509 19.68865 15.45004 13.93727 
##       17       18       19       20 
## 16.56616 15.14733 13.63613 16.26136

Per calcolare i residui invece di usa la funzione resid() passando anche in questo caso lo stesso argomento:

valoriResiduiMult <- resid(linearModelMultiplo)
valoriResiduiMult

##            1            2            3            4            5 
##  0.013546752  0.016115433  0.003765175 -0.004654165 -0.104954889 
##            6            7            8            9           10 
##  0.111941089 -0.003044158  0.088850530 -0.002462277 -0.007942810 
##           11           12           13           14           15 
## -0.004769634 -0.034363753  0.014910594  0.011345859  0.049964672 
##           16           17           18           19           20 
##  0.062733000 -0.066162641 -0.047326362 -0.036130181 -0.061362234

I residui possono essere ottenuti anche attraverso l’attributo “residuals” di linear model:

linearModelMultiplo$residuals

##            1            2            3            4            5 
##  0.013546752  0.016115433  0.003765175 -0.004654165 -0.104954889 
##            6            7            8            9           10 
##  0.111941089 -0.003044158  0.088850530 -0.002462277 -0.007942810 
##           11           12           13           14           15 
## -0.004769634 -0.034363753  0.014910594  0.011345859  0.049964672 
##           16           17           18           19           20 
##  0.062733000 -0.066162641 -0.047326362 -0.036130181 -0.061362234

Dei residui possiamo calcolare mediana, varianza e deviazione standard, mentre non è possibile calcolare il coefficiente di variazione in quanto la media dei residui è 0.
Vediamo dunque nel nostro caso:

median(linearModel$residuals)

## [1] -0.04455862

var(linearModel$residuals)

## [1] 2.906521

sd(linearModel$residuals)

## [1] 1.704852

Anche nel caso multivariato è interessante calcolare i residui standardizzati:

residuiStandardizzatiMulti <- valoriResiduiMult/sd(valoriResiduiMult)
residuiStandardizzatiMulti

##           1           2           3           4           5           6 
##  0.25976522  0.30902087  0.07219897 -0.08924576 -2.01255848  2.14652209 
##           7           8           9          10          11          12 
## -0.05837314  1.70374995 -0.04721531 -0.15230705 -0.09145993 -0.65894083 
##          13          14          15          16          17          18 
##  0.28591753  0.21756209  0.95809566  1.20293426 -1.26869920 -0.90750485 
##          19          20 
## -0.69281290 -1.17664918

I punti ci dicono dove si trovano i residui standardizzati rispetto ai valori stimati con la retta di regressione. In questo caso non si evidenzia nessuna tendenza particolare nel grafico.

Coefficiente di determinazione

Abbiamo già visto in precedenza che il coefficiente di determinazione è definito come il rapporto tra la varianza dei valori stimati con la retta di regressione e la varianza dei valori osservati della variabile dipendente.
L’indice \(D^2\) risulta compreso tra 0 e 1, più è vicino a 1 meglio il modello usato riesce a spiegare i dati.

Per calcolarlo analogamente a prima si può fare:

summary(linearModelMultiplo)$r.square

## [1] 0.9996513

Nel nostro caso quindi, essendo il valore molto vicino ad 1, il modello di regressione spiega in maniera molto significativa i dati.

Screeplot

Lo screeplot è un grafico che rende più semplice la scelta su come partizionare un dendrogramma che si sta analizzando. Sull’asse delle ordinate sono posti il numero di cluster ottenibili e sull’asse delle ascisse le distanze a cui avvengono le aggregazioni.
Se nel passaggio da k a k-1 gruppi la distanza viene incrementata di molto, allora è consigliabile tagliare il dendrogramma in k gruppi.
Si tratta di un metodo empirico ed è consigliabile utilizzare le misure di non omogeneità statistiche (potrebbe suggerire una suddivisione errata).

Costruire uno screeplot è sconsigliato con il metodo del centroide e della mediana (vedremo in seguito che usano i quadrati della distanza) in quanto le agglomerazioni potrebbero verificarsi a livelli di distanza minore o uguale alle precedenti.
Sono invece di aiuto col metodo del legame singolo, completo e medio in cui si utilizza una funzione di distanza.

Analisi del dendrogramma

R mette a disposizione varie funzioni per arricchire il dendrogramma e darci una visione più chiara o per aiutarci ad analizzarlo, vediamole:

• La funzione rect.hclust() ci dà la possibilità di disegnare dei rettangoli intorno ai cluster individuati in base all’altezza a cui vogliamo operare il taglio, oppure specificando il numero di cluster che vogliamo ottenere.

• La funzione cutree() ci consente di ottenere una suddivisione degli individui in cluster in corrispondenza di un livello o di un numero di cluster indicato. Ci restituisce o un vettore o una matrice in cui vengono specificati gli individui in quali cluster sono inseriti.

• La funzione aggregate() ci consente di ricavare misure di sintesi (media, varianza, deviazione standard, …) sui singoli cluster che si sono ottenuti

Ci sono poi le misure più importanti e su cui valuteremo effettivamente la bontà del nostro clustering: misure di non omogeneità statistiche.
Come già detto infatti la traccia di T è univocamente determinata per ogni matrice che descrive p caratteristiche di n individui, allora fissato un numero m di suddivisioni, i cluster devono essere individuati in modo da minimizzare la misura di non omogeneità statistica interna ai cluster, e massimizzare la misura di non omogeneità statistica tra i gruppi.
Considereremo vari metodi gerarchici e vedremo se con lo stesso numero di cluster ci conducono a due diverse partizioni: se capita bisogna scegliere quella che presenta la misura di non omogeneità statistica all’interno dei cluster più piccola (\(trS\)), cioè puntare ad avere una maggiore omogeneità interna ai cluster.

Dato che nel paragrafo successivo si dovranno calcolare le tracce della matrice di non omogeneità statistica tra i cluster, calcoliamo la misura di non omogeneità totale:

numeroRighe <- nrow(matricePopolazionePerTitoloDiStudio)
trHI <- (numeroRighe-1) *sum(
  apply(matricePopolazionePerTitoloDiStudio,2, var ))
trHI

## [1] 1230.762

Vediamo dunque nel dettaglio i vari metodi gerarchici che abbiamo presentato nel paragrafo precedente.

Metodo legame singolo

Il metodo del legame singolo, detto anche nearest neighbour method, individua la distanza tra due cluster come la distanza minima calcolata tra tutte le coppie di individui in cui il primo individuo appartiene al primo cluster, mentre il secondo all’altro cluster preso in considerazione.

Al livello 0 l’algoritmo considera n cluster, uno per ogni individuo.
Al passo 1 si cerca la coppia di individui con la distanza minore e si uniscono in un unico cluster. Si modifica poi la matrice delle distanze scegliendo la distanza come la minima tra quella del primo individuo e quella del secondo individuo del nuovo cluster.
Ad ogni passo dopo che due cluster generici \(G_u\) e \(G_v\) sono stati uniti scegliendo la coppia di cluster meno distante, la distanza tra il nuovo cluster denotato \(G_{uv}\) e un altro cluster \(G_z\) è definita scegliendo dalla precedente matrice delle distanze: \[d_{(uv),z} = \min (d_{uz},d_{vz})\]

Questo metodo ha il vantaggio di essere applicabile a gruppi di qualsiasi forma e di evidenziare la presenza di eventuali valori anomali meglio di altre tecniche, ma ha anche il difetto di basarsi su un singolo legame e non è raro che si possano trovare nello stesso cluster individui piuttosto dissimili: si potrebbero originare delle catene.
Può capitare che due gruppi ben delineati e distinti vengono inseriti nello stesso gruppo erroneamente, dunque non è sempre affidabile il legame singolo.

Effettuiamo il calcolo dei cluster usando il metodo del legame singolo, e vediamo il campo merge dell’output che indica come sono state effettuate le aggregazioni nel corso della computazione:

legameSingolo <- hclust(matriceDistanzeEuclidea,method="single")
legameSingolo$merge

##       [,1] [,2]
##  [1,]   -3  -11
##  [2,]   -1   -2
##  [3,]  -16  -19
##  [4,]   -7   -8
##  [5,]   -9    1
##  [6,]    4    5
##  [7,]   -4   -6
##  [8,]  -15  -18
##  [9,]  -13    6
## [10,]  -10    9
## [11,]    7   10
## [12,]    2   11
## [13,]    3    8
## [14,]  -17   12
## [15,]  -12   14
## [16,]  -14   15
## [17,]   13   16
## [18,]  -20   17
## [19,]   -5   18

Di seguito il dendrogramma e lo screeplot ottenuti utilizzando il metodo:

Valutiamo dunque partizionando in 3, 4 o 5 gruppi come cambia la misura di non omogeneità statistica tra i cluster.

Vediamo come calcoliamo la misura con un partizionamento in 3 gruppi:

taglio <-cutree (legameSingolo , k =3, h = NULL)
num <-table (taglio )
tagliolist <- list(taglio)
agvar <- aggregate (matricePopolazionePerTitoloDiStudio, tagliolist , var)[, -1]

trH1 <-(num [[1]] -1) * sum(agvar [1, ])
if(is.na(trH1))
trH1 <- 0

trH2 <-(num [[2]] -1) *sum(agvar [2, ])
if(is.na(trH2))
trH2 <- 0

trH3 <-(num [[3]] -1) *sum(agvar [3, ])
if(is.na(trH3))
trH3 <- 0

sum <- trH1 + trH2 +trH3
trB <- trHI - sum
trB/trHI

## [1] 0.3212911

Il risultato ci dice che una suddivisione del genere non va affatto bene.

Con 4 gruppi il risultato è decisamente migliore, ci avviciniamo a una suddivisione buona.

## [1] 0.6700402

Con la suddivisione in 5 gruppi non cambia molto rispetto a quanto visto nel caso con 4, comunque sappiamo che il legame singolo si fa influenzare da valori anomali e sappiamo che ce ne sono ben due nel nostro dataset.

## [1] 0.6945096

Metodo legame completo

Il metodo del legame completo, detto anche furthest neighbour method, individua la distanza tra due cluster come la distanza massima calcolata tra tutte le coppie di individui in cui il primo individuo appartiene al primo cluster, mentre il secondo all’altro cluster preso in considerazione.

Al livello 0 l’algoritmo considera n cluster, uno per ogni individuo.
Al passo 1 si cerca la coppia di individui con la distanza minore e si uniscono in un unico cluster. Si modifica poi la matrice delle distanze scegliendo la distanza con gli altri cluster individuata come la maggiore tra quella del primo individuo e quella del secondo individuo del nuovo cluster.
Ad ogni passo dopo che due cluster generici \(G_u\) e \(G_v\) sono stati uniti scegliendo la coppia di cluster meno distante, la distanza tra il nuovo cluster denotato \(G_{uv}\) e un altro cluster \(G_z\) è definita scegliendo dalla precedente matrice delle distanze: \[d_{(uv),z} = \max (d_{uz},d_{vz})\]

Questo metodo è adatto per gruppi che si addensono intorno a un elemento centrale.
Viene privilegiata l’omogeneità dei gruppi e si evita l’effetto catena.
Si nota inoltre che il dendrogramma costruito con questo metodo ha rami più lunghi poiché le distanze sono maggiori. In linea generale il metodo del legame completo è ottimo da utilizzare.

Effettuiamo il calcolo dei cluster usando il metodo del legame completo, e vediamo il campo merge dell’output che indica come sono state effettuate le aggregazioni nel corso della computazione:

legameCompleto <- hclust(matriceDistanzeEuclidea,method="complete")
legameCompleto$merge

##       [,1] [,2]
##  [1,]   -3  -11
##  [2,]   -1   -2
##  [3,]  -16  -19
##  [4,]   -7   -8
##  [5,]   -9    1
##  [6,]   -4   -6
##  [7,]  -15  -18
##  [8,]  -10  -13
##  [9,]  -14  -17
## [10,]    5    8
## [11,]    4    6
## [12,]  -20    3
## [13,]    7    9
## [14,]  -12   10
## [15,]    2   11
## [16,]   13   15
## [17,]   14   16
## [18,]   12   17
## [19,]   -5   18

Di seguito il dendrogramma e lo screeplot ottenuti utilizzando il metodo:

Valutiamo dunque partizionando in 3, 4 o 5 gruppi come cambia la misura di non omogeneità statistica tra i cluster.

Una divisione in 3 gruppi va scartata in quanto non raggiungiamo una misura adeguata.

## [1] 0.5238758

Anche partizionando in 4 cluster non raggiungiamo un risultato positivo, va scartata anche questa suddivisione.

## [1] 0.6390021

Con 5 cluster invece raggiungiamo un risultato ottimo, il migliore ottenuto fino ad ora. Anche dal dendrogramma è evidente che questa suddivisione è ottima.

## [1] 0.8118917

Metodo legame medio

Il metodo del legame medio, detto anche average linkage method, individua la distanza tra due cluster come la media aritmetica delle distanze tra tutte le coppie di individui che compongono due gruppi.

Al livello 0 l’algoritmo considera n cluster, uno per ogni individuo.
Al passo 1 si cerca la coppia di individui con la distanza minore e si uniscono in un unico cluster. Si modifica poi la matrice delle distanze scegliendo la distanza con gli altri cluster: viene individuata come la media calcolata tra i nuovi elementi del cluster (unione degli elementi dei due cluster aggregati) con gli altri insiemi già presenti.
Ad ogni passo dopo che due cluster generici \(G_u\) e \(G_v\) sono stati uniti scegliendo la coppia di cluster meno distante, la distanza tra il nuovo cluster denotato \(G_{uv}\) e un altro cluster \(G_z\) è definita effettuando il calcolo di seguito in base alla matrice delle distanze precedente: \[d_{(uv),z} = \frac{N_u}{N_u + N_v} d_{uz} + \frac{N_v}{N_u + N_v}d_{vz}\]

Dove \(N_u\) e \(N_v\) sono individui nel cluster \(G_u\) e \(G_v\). \(d_{(uv),z}\) rappresenta la misura di distanza media tra gli elementi dei cluster \(G_uv\) e \(G_z\). La procedura si ripete fino ad ottenere un unico cluster formato da tutti gli individui.

Questo metodo ha il seguente svantaggio: se il numero di elementi dei due cluster che si uniscono è molto diverso, la distanza sarà più vicina a quella del cluster più numeroso.

Effettuiamo il calcolo dei cluster usando il metodo del legame medio, e vediamo il campo merge dell’output che indica come sono state effettuate le aggregazioni nel corso della computazione:

legameMedio <- hclust(matriceDistanzeEuclidea,method="average")
legameMedio$merge

##       [,1] [,2]
##  [1,]   -3  -11
##  [2,]   -1   -2
##  [3,]  -16  -19
##  [4,]   -7   -8
##  [5,]   -9    1
##  [6,]   -4   -6
##  [7,]  -15  -18
##  [8,]  -10  -13
##  [9,]    4    5
## [10,]  -14  -17
## [11,]    2    6
## [12,]    3    7
## [13,]    8    9
## [14,]   10   13
## [15,]   11   14
## [16,]  -20   12
## [17,]  -12   15
## [18,]   16   17
## [19,]   -5   18

Di seguito il dendrogramma e lo screeplot ottenuti utilizzando il metodo:

Valutiamo dunque partizionando in 3, 4 o 5 gruppi come cambia la misura di non omogeneità statistica tra i cluster.

Una divisione in 3 gruppi rispetto ai casi precedenti da un risultato migliore, vediamo gli altri.

## [1] 0.6255962

La divisione in 4 gruppi non migliora sostanzialmente la divisione con 3, questo ci fa pensare che il clustering non sia stato effettuato in maniera ottimale.

## [1] 0.684957

Come predetto il risultato con 5 gruppi è deludente, il metodo del legame completo ci ha dato un risultato sicuramente più soddisfacente. in maniera ottimale.

## [1] 0.729401

Metodo del centroide

Il metodo del centroide individua la distanza tra due gruppi come la distanza tra i centroidi, la distanza tra le medie campionarie calcolate sugli individui appartenenti ai due gruppi.

Per questo metodo viene usata la matrice che contiene i quadrati delle singole distanze euclidee.

Questo metodo può portare gruppi di grandi dimensioni a portare dentro di se piccoli gruppi (fenomeni gravitazionali).
Analogamente a prima se uno dei due gruppi uniti ha una numerosità maggiore all’altro, allora il centroide risultante sarà molto vicino a quello del cluster più numeroso.

Calcoliamo dunque la matrice dei quadrati delle singole distanze euclidee che ci servirà per applicare il metodo:

matriceDistanzeEuclideaQuadrato <- matriceDistanzeEuclidea^2

Effettuiamo il calcolo dei cluster usando il metodo del centroide, e vediamo il campo merge dell’output che indica come sono state effettuate le aggregazioni nel corso della computazione:

metodoCentroide <- hclust(matriceDistanzeEuclideaQuadrato,method="centroid")
metodoCentroide$merge

##       [,1] [,2]
##  [1,]   -3  -11
##  [2,]   -1   -2
##  [3,]  -16  -19
##  [4,]   -7   -8
##  [5,]   -9    1
##  [6,]   -4   -6
##  [7,]  -15  -18
##  [8,]  -10  -13
##  [9,]    5    8
## [10,]    4    6
## [11,]  -14  -17
## [12,]    9   11
## [13,]    3    7
## [14,]    2   10
## [15,]   12   14
## [16,]  -20   13
## [17,]  -12   15
## [18,]   16   17
## [19,]   -5   18

In questo caso calcoliamo solo il dendrogramma visto che lo screeplot non è molto utile per questo metodo (e il successivo che vedremo).
Di seguito il dendrogramma dunque:

Valutiamo dunque partizionando in 3, 4 o 5 gruppi come cambia la misura di non omogeneità statistica tra i cluster.

Una divisione in 3 gruppi non dà un risultato molto buono:

## [1] 0.6255962

Anche il risultato con 4 gruppi non è entusiasmante:

## [1] 0.684957

Con 5 gruppi raggiungiamo una soglia buona, ma col metodo del legame completo abbiamo raggiunto un risultato migliore:

## [1] 0.729401

Metodo della mediana

Infine è stato utilizzato anche il metodo della mediana senza nessun risultato significativo.

metodoMediana <- hclust(matriceDistanzeEuclideaQuadrato,method="median")
metodoMediana$merge

##       [,1] [,2]
##  [1,]   -3  -11
##  [2,]   -1   -2
##  [3,]  -16  -19
##  [4,]   -7   -8
##  [5,]   -9    1
##  [6,]   -4   -6
##  [7,]  -15  -18
##  [8,]  -10  -13
##  [9,]    5    8
## [10,]    4    6
## [11,]  -14  -17
## [12,]    9   11
## [13,]    3    7
## [14,]    2   10
## [15,]   12   14
## [16,]  -20   13
## [17,]  -12   15
## [18,]   16   17
## [19,]   -5   18

Il corrispondente dendrogramma è il seguente:

Di seguito la misura di non omogeneità statistica ottenuta con 3, 4 e 5 cluster:

## [1] 0.6255962

## [1] 0.684957

## [1] 0.729401

In conclusione da questa analisi abbiamo individuato che il clustering migliore lo otteniamo utilizzando il metodo del legame completo con 5 cluster.
Il clustering ottenuto è dunque:

##              Piemonte         Valle d'Aosta               Liguria 
##                     1                     1                     2 
##             Lombardia   Trentino Alto Adige                Veneto 
##                     1                     3                     1 
## Friuli-Venezia Giulia        Emilia-Romagna               Toscana 
##                     1                     1                     2 
##                Umbria                Marche                 Lazio 
##                     2                     2                     2 
##               Abruzzo                Molise              Campania 
##                     2                     4                     4 
##                Puglia            Basilicata              Calabria 
##                     5                     4                     4 
##               Sicilia              Sardegna 
##                     5                     5

\(C_1 = \{Piemonte, Valle D'aosta, Lombardia, Veneto, Friuli-Venezia\ Giulia,Emilia-Romagna\}\)
\(C_2 = \{Liguria,Toscana,Umbria,Marche,Lazio,Abruzzo\}\)
\(C_3 = \{Trentino Alto Adige\}\)
\(C_4 = \{Molise,Campania,Basilicata,Calabria\}\)
\(C_5 = \{Puglia,Sicilia,Sardegna\}\)

La misura di non omogeneità statistica tra cluster ottenuta è pari a:

## [1] 0.8118917

Passiamo dunque all’analisi con i metodi non gerarchici per migliorare ulteriormente il partizionamento ottenuto.