Media, mediana e deviazione standard in R

Media e mediana

I due riepiloghi più usati di tutta la statistica richiedono una chiamata di funzione ciascuno:

mpg <- mtcars$mpg

print(mean(mpg))
print(median(mpg))

mean() somma tutto e divide per il numero di valori. median() ordina i valori e prende quello centrale (o fa la media dei due centrali quando il numero è pari). Per le 32 auto di mtcars, il consumo medio è di circa 20.1 mpg e la mediana è 19.2: valori vicini, il che ti dice che i dati non sono molto asimmetrici. Quando invece differiscono molto, anche questa è un'informazione; ci arriveremo alla fine.

Una cosa che sorprende tutti: se il vettore contiene anche un solo NA, entrambe le funzioni restituiscono NA:

x <- c(4, 8, NA, 6)

print(mean(x))
print(mean(x, na.rm = TRUE))

È voluto: R si rifiuta di fingere in silenzio che il valore mancante non ci sia. na.rm = TRUE significa "calcola sui valori che hai". Quasi tutte le funzioni di riepilogo di questo articolo lo accettano. La storia completa di come si propaga NA è nei valori mancanti.

Deviazione standard e varianza

sd() misura la dispersione: più o meno, quanto un valore tipico si allontana dalla media, nelle stesse unità dei dati. var() è il suo quadrato:

mpg <- mtcars$mpg

print(sd(mpg))
print(var(mpg))
print(sd(mpg)^2)

Una deviazione standard di circa 6 mpg significa che le auto stanno tipicamente entro circa 6 mpg dalla media di 20.1. Siccome sd() è nelle unità dei dati, è quella che riporti; var() compare soprattutto dentro altre formule.

Ecco il dettaglio che conta per gli esami: sd() e var() calcolano la statistica campionaria, cioè dividono la somma dei quadrati degli scarti per n − 1, non per n:

x <- c(2, 4, 4, 4, 5, 5, 7, 9)
n <- length(x)

sample_var <- sum((x - mean(x))^2) / (n - 1)
population_var <- sum((x - mean(x))^2) / n

print(sample_var)      # what var(x) gives you
print(var(x))
print(population_var)  # divide by n instead

Perché n − 1? Perché hai stimato la media dagli stessi dati, e gli scarti attorno a quella media stimata sono sistematicamente un po' troppo piccoli; dividere per n − 1 corregge questo effetto. Siccome i tuoi dati sono quasi sempre un campione di qualcosa di più grande, la versione con n − 1 è quella che ti serve. Se hai davvero l'intera popolazione (ogni studente della classe, ogni prodotto del catalogo), moltiplica: var(x) * (n - 1) / n.

Errore standard della media

Deviazione standard ed errore standard si confondono di continuo, quindi tienili separati: la sd descrive i dati, l'errore standard descrive la tua stima della media. R non ha una se() integrata, ma la formula è una riga:

mpg <- mtcars$mpg

se <- sd(mpg) / sqrt(length(mpg))
print(se)

L'errore standard si riduce quando il campione cresce (raccogli quattro volte i dati e si dimezza), perché un campione più grande individua la media con più precisione. La sd invece non si riduce con la dimensione del campione: le auto sono varie quanto sono, indipendentemente da quante ne misuri. L'errore standard è il mattone degli intervalli di confidenza, ed è lì che si guadagna il posto.

summary(): la panoramica in una chiamata

summary() ti dà il riepilogo a cinque numeri più la media in un colpo solo, e funziona su interi data frame:

print(summary(mtcars$mpg))

print(summary(mtcars[, c("mpg", "wt", "hp")]))

Chiamata su un data frame, riassume ogni colonna: le colonne numeriche ottengono minimo, quartili, media e massimo, e i factor i conteggi per livello. È la prima cosa da eseguire su qualsiasi dataset appena caricato: i valori impossibili (un'età negativa, un massimo di 9999) saltano subito all'occhio.

Quantili, intervallo e IQR

quantile() generalizza la mediana a qualsiasi punto di taglio:

mpg <- mtcars$mpg

print(quantile(mpg))
print(quantile(mpg, probs = c(0.1, 0.5, 0.9)))

print(range(mpg))
print(IQR(mpg))
print(max(mpg) - min(mpg))

Senza argomenti restituisce minimo, quartili e massimo. Passa probs = per punti di taglio specifici: il 10° e il 90° percentile qui sopra delimitano la zona in cui vive la maggior parte dei dati. IQR() (la distanza tra il 25° e il 75° percentile) è una misura di dispersione che, a differenza di sd(), non si lascia trascinare dagli outlier. range() restituisce minimo e massimo come coppia.

La moda: mode() di R NON fa questo

Questa frega tutti esattamente una volta. R ha una funzione chiamata mode(), e non ha niente a che fare con la statistica: riporta il tipo di memorizzazione di un oggetto:

x <- c(3, 7, 7, 2, 7, 5, 3)

print(mode(x))   # "numeric" - not what you wanted

counts <- table(x)
print(counts)
print(names(which.max(counts)))

L'idioma da ricordare: table(x) conta quante volte compare ogni valore, which.max() trova il conteggio più grande e names() estrae il valore stesso. Nota che torna come stringa di caratteri (i nomi di una table lo sono sempre); avvolgilo in as.numeric() se devi usarlo nei calcoli. Se due valori sono a pari merito, which.max() restituisce in silenzio solo il primo: controlla tu la tabella quando i pareggi sono plausibili.

Media o mediana: quale riportare

La media usa ogni valore, ed è la sua forza e la sua debolezza: un solo valore estremo la trascina. Alla mediana interessa solo il centro, quindi gli outlier la toccano appena:

salaries <- c(42000, 45000, 47000, 48000, 51000, 54000, 58000)
print(mean(salaries))
print(median(salaries))

# The CEO joins the dataset:
salaries <- c(salaries, 900000)
print(mean(salaries))
print(median(salaries))

Un solo valore aggiunto spinge la media da circa 49.300 a oltre 155.000, un numero che non descrive nessuno nei dati, mentre la mediana si sposta solo da 48.000 a 49.500. Ecco perché redditi, prezzi delle case e durate dei ricoveri si riportano come mediane: con dati asimmetrici e una coda lunga la media diventa fuorviante. Per dati più o meno simmetrici le due coincidono e la media va benissimo (ed è statisticamente più efficiente). Un rapido istogramma ti dice in quale situazione ti trovi, e confrontare la media con la mediana è di per sé un controllo dell'asimmetria in una riga.

Cosa ti porti a casa

• mean(x) e median(x): aggiungi na.rm = TRUE quando ci sono valori mancanti.
• sd(x) e var(x) calcolano la statistica campionaria (il denominatore n − 1), ed è ciò che ti serve.
• L'errore standard è sd(x) / sqrt(length(x)): misura quanto bene conosci la media, non quanto sono dispersi i dati.
• summary() su un data frame appena caricato è il controllo di buon senso più veloce di R.
• La moda è names(which.max(table(x))); mode() di R riguarda i tipi di memorizzazione.
• Dati asimmetrici o outlier: riporta la mediana. Dati simmetrici: la media va bene.

Prossimo passo: misurare come due variabili si muovono insieme, con la correlazione tramite cor() e cor.test().