6.1 Alcune funzioni di R per la statistica descrittiva
6.1.1 Analisi descrittive univariate
6.1.1.1 Una singola variabile categorica
Nel caso ci interessi sintetizzare le informazioni relative ad una variabile categorica, gli strumenti che possiamo utilizzare sono la distribuzione di frequenze insieme ad una rappresentazione grafica della stessa.
Prima di mostrare alcuni esempi, vi ricordiamo che in R le variabili categoriche corrispondono a vettori di tipo factor, che abbiamo presentato nella Sezione 1.3.4. Se la variabile che intendete analizzare non fosse codificata come factor, dovete prima convertirla in quel formato usando la funzione as.factor(). Ad esempio, per convertire la variabile numerica IndustryCode nel data frame forbes94 dobbiamo eseguire il codice seguente:
> load("Forbes94.RData")
> forbes94$IndustryCode <- as.factor(forbes94$IndustryCode)
> str(forbes94$IndustryCode) Factor w/ 20 levels "0","1","2","3",..: 12 12 12 12 12 12 12 12 12 12 ...Ricordate che l’operatore $ serve per riferirsi ad una colonna di un data frame.
Per costruire la distribuzione di frequenze di una variabile categorica in R dobbiamo usare la funzione table(), già descritta nella Sezione 2.1. Ad esempio, la distribuzione di frequenze della variabile categorica WideIndustry si ottiene con il codice seguente:
> table(forbes94$WideIndustry)
Aerospacedefense Business Capital goods
19 27 19
Chemicals ComputersComm Construction
25 67 11
Consumer Energy Entertainment
54 42 27
Financial Food Forest
168 62 19
Health Insurance Metals
49 54 18
Retailing Transport Travel
46 15 15
Utility
63 Come sempre, è possibile memorizzare il risultato della funzione table() in un oggetto (si tratterebbe di un semplice vettore numerico contenente le frequenze assolute e i cui elementi sarebbero chiamati come i livelli, ovvero le modalità, della variabile) e utilizzarlo successivamente per ulteriori analisi.
Per produrre la distribuzione di frequenze relative possiamo usare la funzione prop.table(), la quale richiede in input la tabella con le frequenze assolute:
> absfreq <- table(forbes94$WideIndustry)
> prop.table(absfreq)
Aerospacedefense Business Capital goods
0.02375 0.03375 0.02375
Chemicals ComputersComm Construction
0.03125 0.08375 0.01375
Consumer Energy Entertainment
0.06750 0.05250 0.03375
Financial Food Forest
0.21000 0.07750 0.02375
Health Insurance Metals
0.06125 0.06750 0.02250
Retailing Transport Travel
0.05750 0.01875 0.01875
Utility
0.07875 Una volta creata la distribuzione di frequenze, possiamo rappresentarla graficamente con un diagramma a torta o a barre. In R questi si ottengono rispettivamente con le funzioni pie() (già descritta nella Sezione 2.1) e barplot(), le quali richiedono in input la distribuzione di frequenze (assolute o relative):
> pie(absfreq, main = "Distribuzione di WideIndustry") # diagramma a torta
> barplot(absfreq, main = "Distribuzione di WideIndustry") # diagramma a barreFigura 6.1: Diagramma a torta e a barre con R.
6.1.1.2 Una singola variabile numerica
6.1.1.2.1 Distribuzione di frequenze
La distribuzione di una variabile numerica è rappresentata da una tabella che riporta le frequenze con cui sono stati osservati i diversi valori, nel caso la variabile sia discreta, o le classi in cui è stata divisa, nel caso in cui la variabile sia continua. In entrambi i casi è possibile utilizzare ancora la funzione table(), ma nel caso di variabili continue sarà prima necessario dividere i valori osservati in classi. Una funzione utile per generare tali classi è cut(), che genera un vettore factor con le classi specificate e che include i seguenti argomenti:
x, il vettore numerico con i valori di cui vogliamo costruire le classibreaks, che consente di specificare quali classi intendiamo usare; può essere specificato nei seguenti modi alternativi- un singolo valore numerico che indica il numero di classi di uguale ampiezza che vogliamo usare
- un vettore di valori numerici che identificano gli estremi delle classi (eventualmente di diversa ampiezza) da usare
right, un valore logico che indica se le classi devono essere chiuse a sinistra (FALSE) o a destra (TRUE)labels, un vettore di etichette da utilizzare per identificare le classi; se fissato pari aFALSE, le classi saranno identificate da dei numeri interi
Per esempio, supponiamo di voler costruire la distribuzione di frequenze della variabile Salary, che fornisce lo stipendio pagato ai CEO inclusi nel data set nel 1994, usando 5 classi di uguale ampiezza chiuse a sinistra. Il seguente codice permette di creare una nuova variabile, che decidiamo di chiamare Salary_c e che aggiungiamo al data frame esistente:
> forbes94$Salary_c <- cut(forbes94$Salary, breaks = 5, right = FALSE)
> table(forbes94$Salary_c)
[1.59e+04,5.55e+05) [5.55e+05,1.09e+06)
382 376
[1.09e+06,1.63e+06) [1.63e+06,2.16e+06)
24 4
[2.16e+06,2.7e+06)
3 da cui si vede chiaramente che la distribuzione dei salari per questi CEO è decisamente asimmetrica a destra.
Se invece volessimo utilizzare le classi \([0, 200000)\), \([200000, 500000)\), \([500000, 800000)\), \([800000, 1500000)\) e \([1500000, 3000000]\), allora dovremmo utilizzare la funzione cut() nel modo seguente:
> forbes94$Salary_c <- cut(forbes94$Salary,
+ breaks = c(0, 200000, 500000, 800000, 1500000, 3000000), right = FALSE)
> table(forbes94$Salary_c)
[0,2e+05) [2e+05,5e+05) [5e+05,8e+05)
13 268 351
[8e+05,1.5e+06) [1.5e+06,3e+06)
147 10 Per rappresentare graficamente la distribuzione di una variabile numerica sono disponibili vari tipi di grafici tra cui il diagramma ad aste, l’istogramma e il box-plot.
6.1.1.2.2 Diagramma ad aste
Possiamo creare un diagramma ad aste usando la funzione plot() direttamente sul risultato fornito da table() e fissando l’argomento type al valore "h". Ad esempio, il diagramma ad aste per la variabile Age è ottenuto come segue:
> plot(table(forbes94$Age), type = "h")
Figura 6.2: Diagramma ad aste per la variabile Age.
L’aspetto finale del grafico può essere modificato attraverso gli argomenti della funzione plot(). Se ad esempio volessimo aggiungere le etichette agli assi dovremmo usare gli argomenti xlab e ylab, ovvero
> plot(table(forbes94$Age), type = "h",
+ xlab = "Age", ylab = "Frequenza assoluta",
+ main = "Diagramma ad aste della variabile Age")
Figura 6.3: Diagramma ad aste per la variabile Age.
6.1.1.2.3 Istogramma
Gli istogrammi sono rappresentazioni grafiche della distribuzione di frequenze di una variabile divisa per classi. In R è possibile creare un istogramma con la funzione hist(). Come è noto, è possibile creare un istogramma usando classi di uguale o diversa ampiezza. La funzione hist() permette di specificare tale scelta attraverso l’argomento breaks, il cui funzionamento è identico al medesimo argomento descritto per la funzione cut(), ovvero:
- se non specifichiamo il valore di
breaks, R utilizza alcune regole interne per sceglierne un valore ragionevole; consigliamo di usare questa opzione, perché funziona bene nella maggioranza dei casi, a meno che non ci siano motivi particolari per usare una regola diversa breakspuò essere specificato come singolo valore numerico (intero positivo), il quale indicherà il numero (indicativo45) di classi di uguale ampiezza che vogliamo usarebreakspuò anche essere specificato come vettore numerico, i cui elementi corrisponderanno agli estremi delle classi che intendiamo usare
In aggiunta a breaks, la funzione hist() contiene molti altri argomenti con cui è possibile modificare l’aspetto finale del grafico, tra cui main, xlab e ylab, che funzionano in modo analogo a quanto visto finora.
Vediamo ora alcuni esempi di istogrammi per la variabile Salary:
> hist(forbes94$Salary, main = "Istogramma di Salary", xlab = "Salary (in $)",
+ ylab = "Frequenza assoluta")
Figura 6.4: Esempio di istogramma per la variabile Salary (classi di uguale ampiezza).
> hist(forbes94$Salary, main = "Istogramma di Salary", xlab = "Salary (in $)",
+ ylab = "Frequenza assoluta", breaks = 50)
Figura 6.5: Esempio di istogramma per la variabile Salary (classi di uguale ampiezza).
> hist(forbes94$Salary, main = "Istogramma di Salary", xlab = "Salary (in $)",
+ ylab = "Densità di frequenze",
+ breaks = c(0, 100000, 300000, 350000, 500000, 800000, 1000000, 3000000))
Figura 6.6: Esempio di istogramma per la variabile Salary (classi di diversa ampiezza).
Questi esempi mostrano chiaramente che la distribuzione di Salary è asimmetrica a destra. Vi facciamo inoltre notare che nell’ultimo esempio, in cui abbiamo utilizzato 7 classi di diversa ampiezza, R ha automaticamente (e correttamente!) deciso di riportare sull’asse verticale la densità di frequenze, mentre nei due esempi precedenti sull’asse verticale ha riportato le frequenze assolute delle classi. Anche nel caso di classi di uguale ampiezza è possibile costruire l’istogramma con le densità fissando a FALSE il valore dell’argomento freq46.
6.1.1.2.4 Box-plot
In aggiunta al diagramma ad aste e all’istogramma, è utile produrre anche il box-plot, che consente di integrare le informazioni fornite dagli altri grafici. Come abbiamo già visto nella Sezione 2.1, in R è possibile costruire un box-plot attraverso la funzione boxplot(). Ricordiamo che il box-plot della variabile Salary si ottiene come segue:
> boxplot(forbes94$Salary)Figura 6.7: Esempio di box-plot con R.
Un utile commento finale su boxplot() riguarda la possibilità di salvare output aggiuntivo rispetto al grafico. In particolare, boxplot() consente di salvare in una lista (vedi la Sezione 1.3.6) i valori degli indici che servono per costruire il grafico. Per esempio, per la variabile Salary otteniamo
> bp_out <- boxplot(forbes94$Salary)
> str(bp_out)List of 6
$ stats: 'integer' num [1:5, 1] 18600 433333 569231 750000 1200000
$ n : num 789
$ conf : num [1:2, 1] 551419 587043
$ out : num [1:18] 1750000 1400000 1500000 1375000 2000000 ...
$ group: num [1:18] 1 1 1 1 1 1 1 1 1 1 ...
$ names: chr "1"Come vedete, bp_out è una lista che contiene 6 elementi (ovvero altri oggetti), tra cui stats, che contiene i valori dei 5 indici necessari per costruire il box-plot (minimo, primo quartile, mediana, terzo quartile e massimo), e out, il vettore delle osservazioni outlier. Concludiamo quindi che nella distribuzione di Salary risultano essere presenti 18 outlier:
[1] 1750000 1400000 1500000 1375000 2000000 1900000
[7] 1250000 1483500 1454000 2633570 1240000 1841000
[13] 1231730 2667730 2700000 1525000 1397810 15049406.1.1.2.5 Indici di sintesi
Le principali funzioni disponibili in R per il calcolo di indici di sintesi per una variabile numerica sono:
min(), per il calcolo del minimo valore osservato di una variabilemax(), per il calcolo del massimo valore osservato di una variabilemean(), per il calcolo della media campionariamedian(), per il calcolo della mediana campionariavar(), per il calcolo della varianza campionariasd(), per il calcolo della deviazione standard campionariaquantile(), per il calcolo dei quantili campionarisummary(), la quale fornisce in automatico una lista con alcuni degli indici sopra citati
Nel caso la variabile contenesse dei valori mancanti (ovvero NA), queste funzioni ad eccezione di summary() richiedono venga fissato a TRUE il valore dell’argomento na.rm, tramite il quale indichiamo ad R di rimuovere internamente (ma non dal data set) i dati mancanti prima di effettuare il calcolo. Vediamo un esempio usando sempre la variabile Salary:
> min(forbes94$Salary, na.rm = TRUE)[1] 18600> max(forbes94$Salary, na.rm = TRUE)[1] 2700000> mean(forbes94$Salary, na.rm = TRUE)[1] 613458> median(forbes94$Salary, na.rm = TRUE)[1] 569231> var(forbes94$Salary, na.rm = TRUE)[1] 7.858e+10> sd(forbes94$Salary, na.rm = TRUE)[1] 280319> quantile(forbes94$Salary, na.rm = TRUE) 0% 25% 50% 75% 100%
18600 433333 569231 750000 2700000 > summary(forbes94$Salary) Min. 1st Qu. Median Mean 3rd Qu. Max.
18600 433333 569231 613458 750000 2700000
NA's
11 La funzione quantile() consente di calcolare i percentili di un set di dati numerici. L’ordine dei percentili che vogliamo calcolare viene specificato attraverso l’argomento probs. Per esempio, i percentili di ordine 10 e 95 per la variabile Salary sono calcolati come segue:
> quantile(forbes94$Salary, na.rm = TRUE, probs = c(0.10, 0.95)) 10% 95%
325149 1025640 6.1.2 Analisi descrittive bivariate
6.1.2.1 Due variabili categoriche
Nel caso volessimo analizzare congiuntamente due variabili categoriche, lo strumento principale da utilizzare consiste nella tabella a doppia entrata. In R possiamo costruire una tabella a doppia entrata con la funzione table(), che abbiamo già incontrato in varie occasioni. Per ottenere una tabella a doppia entrata bisogna però specificare due variabili, rispettivamente la variabile di riga e quella di colonna (in quest’ordine). Il codice seguente permette di costruire la tabella a doppia entrata per la variabile MBA contro WideIndustry:
> table(forbes94$WideIndustry, forbes94$MBA)
0 1
Aerospacedefense 15 4
Business 22 5
Capital goods 13 6
Chemicals 17 8
ComputersComm 46 21
Construction 7 4
Consumer 36 18
Energy 32 10
Entertainment 23 4
Financial 112 56
Food 45 17
Forest 12 7
Health 37 12
Insurance 43 11
Metals 14 4
Retailing 43 3
Transport 11 4
Travel 13 2
Utility 48 15La rappresentazione grafica di una tabella a doppia entrata tramite un diagramma a barre sovrapposte si ottiene con la funzione plot() in cui le variabili devono essere specificate attraverso una formula:
> plot(MBA ~ WideIndustry, data = forbes94)
Figura 6.8: Grafico a barre sovrapposte per la variabile MBA rispetto a WideIndustry.
6.1.2.2 Due variabili numeriche
In R è possibile creare un diagramma di dispersione con la funzione plot() in cui specifichiamo come prima la variabile \(x\) e come seconda la \(y\)47. Come per le altre funzioni grafiche, possiamo aggiungere anche altri argomenti con i quali possiamo modificare l’aspetto finale del grafico. Il codice che segue crea il diagramma di dispersione per la variabile Salary contro la variabile Age:
> plot(forbes94$Age, forbes94$Salary, xlab = "Age (in years)",
+ ylab = "Salary (in $)", main = "Diagramma di dispersione")
Figura 6.9: Diagramma di dispersione per la variabile Salary rispetto a Age.
Per il calcolo della covarianza e dell’indice di correlazione lineare campionari possiamo usare rispettivamente le funzioni cov() e cor(). Nel caso includessimo in queste funzioni più di due variabili, queste restituirebbero in output rispettivamente le matrici delle covarianze e delle correlazioni lineari campionarie. Se alcune delle variabili contengono dei dati mancanti, è opportuno includere in queste funzioni anche l’argomento use fissandolo al valore "complete.obs"48. Per calcolare le covarianze e le correlazioni tra le variabili Age, Salary, Bonus e Other possiamo eseguire il codice che segue:
> cov(forbes94[, c("Age", "Salary", "Bonus", "Other")],
+ use = "complete.obs") Age Salary Bonus Other
Age 41.59 4.059e+05 7.083e+05 4.465e+05
Salary 405931.80 6.886e+10 6.194e+10 9.696e+10
Bonus 708346.61 6.194e+10 8.676e+11 4.023e+11
Other 446454.17 9.696e+10 4.023e+11 1.671e+12> cor(forbes94[, c("Age", "Salary", "Bonus", "Other")],
+ use = "complete.obs") Age Salary Bonus Other
Age 1.00000 0.2399 0.1179 0.05356
Salary 0.23987 1.0000 0.2534 0.28583
Bonus 0.11792 0.2534 1.0000 0.33416
Other 0.05356 0.2858 0.3342 1.00000Le correlazioni tra queste variabili sono tutte positive, ma comunque deboli o molto deboli. Ciò è dovuto principalmente alla presenza di numerosi outlier e all’eterogeneità dei valori osservati.
6.1.3 Ulteriori rappresentazioni grafiche
6.1.3.1 Il diagramma di Pareto
Il diagramma di Pareto è una rappresentazione grafica utile per separare le modalità di una variabile categorica che sono state osservate più frequentemente da quelle meno frequenti. Purtroppo non c’è nessuna funzione nei pacchetti base di R che consente di produrre un diagramma di Pareto. Tuttavia, è possibile scaricare e installare un pacchetto aggiuntivo, chiamato qcc, che contiene alcune funzioni dedicate al controllo statistico della qualità in produzione, tra cui anche una funzione per disegnare il diagramma di Pareto. Possiamo installare il pacchetto seguendo le istruzioni riportate in Sezione 1.5 o in alternativa eseguendo il codice seguente:
> install.packages("qcc") # scarica il pacchetto qcc
> library(qcc) # carica il pacchetto qccIl pacchetto qcc contiene una funzione chiamata pareto.chart() che consente di produrre direttamente il diagramma di Pareto. Questa funzione richiede che venga fornito in input un vettore numerico con le frequenze assolute con cui si sono osservate le varie modalità. Per ottenere un grafico più informativo, consigliamo anche di assegnare dei nomi agli elementi del vettore, i quali saranno riportati sull’asse orizzontale del grafico. Il seguente codice dapprima carica alcuni dati contenuti nel data frame pareto_dat e relativi al numero di pezzi difettosi rilevati in un’azienda (colonna Total), insieme alla corrispondente causa di difettosità (colonna Cause), e successivamente utilizza questi dati per costruire il diagramma di Pareto:
> load("pareto_data.RData")
> pdat <- pareto_dat$Total
> names(pdat) <- pareto_dat$Cause
> pareto.chart(data = pdat, main = "Diagramma di Pareto")
Figura 6.10: Esempio di diagramma di Pareto ottenuto con la funzione pareto.chart() del pacchetto qcc.
Pareto chart analysis for pdat
Frequency Cum.Freq.
Compressor 15.000 15.000
Machine adjustment 11.000 26.000
Packing 10.000 36.000
Supplier 9.000 45.000
Operation conditions 8.000 53.000
Transport agency 3.000 56.000
Measurement apparatus 3.000 59.000
Receptionist 2.000 61.000
Transport method 2.000 63.000
Recording method 2.000 65.000
Recording Operator 1.000 66.000
Storage operators 1.000 67.000
Raw materials reception 1.000 68.000
Pareto chart analysis for pdat
Percentage Cum.Percent.
Compressor 22.059 22.059
Machine adjustment 16.176 38.235
Packing 14.706 52.941
Supplier 13.235 66.176
Operation conditions 11.765 77.941
Transport agency 4.412 82.353
Measurement apparatus 4.412 86.765
Receptionist 2.941 89.706
Transport method 2.941 92.647
Recording method 2.941 95.588
Recording Operator 1.471 97.059
Storage operators 1.471 98.529
Raw materials reception 1.471 100.0006.1.3.2 Rappresentazione grafica di una serie storica
Un’altra situazione che si incontra di frequente quando si analizzano dei dati è quella in cui le osservazioni sono ordinate temporalmente. Questo tipo di dati vengono chiamati serie storiche. La funzione ts() consente di creare una serie storica e di definirne le caratteristiche in termini di numero di osservazioni per anno, data di inizio e data di fine. A titolo di esempio, carichiamo il file advertising_sales.RData, che contiene il data frame advertising_sales con le serie storiche mensili delle vendite e delle spese in pubblicità (entrambe in migliaia di dollari) di un prodotto per il fitness. I dati disponibili sono stati raccolti a partire dal mese di gennaio 2012. Il seguente codice permette di definire i dati come serie storiche applicando la funzione ts() con l’argomento frequency fissato a 12 (ovvero 12 osservazioni in un anno) e l’argomento start al vettore c(2012, 1), che indicano rispettivamente l’anno e il mese a cui fa riferimento la prima osservazione. Infine, con la funzione plot() creiamo il grafico temporale delle due serie, che rappresentiamo con colori e tipo di linee (argomenti col e lty) diversi49:
> load("advertising_sales.RData")
> adv_sales <- ts(advertising_sales[, c(2, 3)], frequency = 12, start = c(2012, 1))
> plot(adv_sales, plot.type = "single", col = c("red", "blue"),
+ lwd = 2, lty = c(1, 2), ylab = "", main = "Grafico per serie storiche",
+ sub = "(linea continua = spese in pubblicità, linea tratteggiata = vendite)")Figura 6.11: Esempio di grafico per serie storiche.
6.1.3.3 Come verificare se una variabile è distribuita in modo normale
Molti modelli statistici richiedono che le popolazioni di riferimento siano distribuite seconda una normale. Pertanto, è necessario verificare che questa ipotesi sia almeno approssimativamente soddisfatta dai dati campionari che stiamo analizzando. In caso contrario, i risultati che otterremo dalle analisi potrebbero essere distorti e non attendibili.
Uno strumento semplice ma efficace per verificare la normalità di una distribuzione di un set di dati numerici è il normal probability plot, che si può creare in R attraverso la funzione qqnorm(). Questa funzione richiede solo l’indicazione della variabile da usare, oltre eventualmente ad altri argomenti per migliorare l’aspetto del grafico (main, xlab, ylab, col, ecc.).
Dopo aver eseguito qqnorm(), la funzione qqline() consente di aggiungere al grafico una retta che aiuta a valutare la vicinanza della distribuzione empirica dei dati alla distribuzione normale. Anche questa funzione richiede obbligatoriamente di indicare solamente il nome della variabile oggetto dell’analisi.
Il seguente codice genera il normal probability plot per la variabile Salary:
> qqnorm(forbes94$Salary, main = "Normal probability plot di Salary")
> qqline(forbes94$Salary)
Figura 6.12: Esempio di normal probability plot per la variabile Salary.
Dal grafico si conclude che la distribuzione di Salary è asimmetrica a destra e pertanto non si può ritenere distribuita come una normale.
Come secondo esempio consideriamo la variabile Age:
> qqnorm(forbes94$Age, main = "Normal probability plot di Age")
> qqline(forbes94$Age)
Figura 6.13: Esempio di normal probability plot per la variabile Age.
In questo caso la distribuzione empirica di Age è molto più vicina, anche se non perfettamente, ad una distribuzione normale50.
Come è indicato anche nell’help della funzione
hist(), infatti, il valore che forniamo per il numero di classi di uguale ampiezza è solo un “suggerimento”, ma internamente tale valore viene rivisto in funzione di alcune regole, che non è nostro obiettivo descrivere ora.↩Se provaste a fissare
freq = TRUEnel caso di classi di diversa ampiezza, otterreste un istruttivo messaggio di errore.↩In alternativa, è possibile specificare le variabili attraverso una formula, ad esempio
plot(Salary $\sim$ Age, data = forbes94).↩Sono disponibili altre scelte per l’argomento
use, ma non le presentiamo in questo manuale. Vi invitiamo comunque a leggere l’help di queste funzioni.↩Gli argomenti
lwd,ltyesubservono per definire lo spessore della linea (2significa linee spesse il doppio del normale), il tipo di linea (1significa una linea continua, mentre 2 una linea tratteggiata) e un sottotitolo, riportato nella parte bassa del grafico.↩Ricordiamo che il normal probability plot è uno strumento grafico e pertanto aiuta a discriminare bene se ci troviamo vicino alla situazione di normalità o meno solo nei casi più evidenti. In molte situazioni il grafico non consente di arrivare ad una conclusione univoca sulla normalità. In questi casi è necessario utilizzare strumenti più elaborati che non sono oggetto di discussione né qui né nel corso.↩