Loop e vettorializzazione con R

Questo blog ha come obiettivo quello di fornire gli strumenti per "rompere il ghiaccio" con R, lasciando poi agli eventuali interessi e all'iniziativa personale l'approfondimento del linguaggio di programmazione, che ciascuno potrà realizzare in base ai proprio bisogni. Tuttavia occasionalmente può essere interessante affrontare qualche argomento un poco più tecnico, come quello che vediamo oggi.

Il punto di partenza è rappresentato dal calcolo delle statistiche elementari parametriche. Per un campione che include n dati:

a) la media, la misura di posizione dei dati, viene calcolata come

b) la deviazione standard, la misura di dispersione dei dati, viene calcolata come

essendo

la varianza dei dati.

Invece di calcolare le statistiche elementari parametriche servendoci di qualcuna delle funzioni disponibili in R le vogliamo calcolare seguendo alla lettera la loro definizione matematica. Come si vede i calcoli sono elementari, e di fatto l'unico problema da risolvere in termini di programmazione è trovare un metodo per il calcolo della sommatoria Σ  che ci serve per ottenere:

1) la somma di tutti i dati

da dividere poi per il numero n dei dati al fine di calcolare la media;

2) la somma dei quadrati delle differenze tra ciascun dato e la media, ovvero la devianza

da dividere poi per il numero dei dati meno uno per calcolare la varianza dalla quale otteniamo infine la deviazione standard.

Il primo metodo per il calcolo di una sommatoria che impieghiamo fa ricorso ad una iterazione (o ciclo o loop), una struttura di controllo che consente di eseguire ripetutamente una o più istruzioni, fino al verificarsi di una condizione specificata.

Accertatevi innanzitutto di avere installato il pacchetto psychTools che contiene i dati che ci servono, o in alternativa procedete come indicato nel post Il set di dati galton nel quale trovate anche illustrato il significato dei dati impiegati. Quindi copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# STATISTICHE ELEMENTARI PARAMETRICHE con loop
#
library(psychTools) # carica il pacchetto che include il set di dati galton

x <- galton$parent # salva nel vettore x le altezze dei padri
#
n <- length(x) # calcola il numero dei dati

#

### (a)

#
sx <- 0 # inizializza la sommatoria dei dati
for (i in 1:n) {
    sx <- sx+x[i] # calcola la sommatoria dei dati
}

#
media <- sx/n # calcola la media

#
### (b)

#
devianza <- 0 # inizializza la sommatoria degli scarti quadratici (devianza)
#
for (i in 1:n) {
    devianza <- devianza+(x[i]-media)^2 # calcola la sommatoria degli scarti quadratici
}
#
varianza <- devianza/(n-1) # calcola la varianza
ds <- sqrt(varianza) # calcola la deviazione standard

#
### (c)

#
statistica <- c("numero dei dati =", "media =", "varianza =", "deviazione standard =") # predispone i nomi delle statistiche calcolate 
risultato <- c(n, media, varianza, ds) # predispone i risultati delle statistiche calcolate
print(data.frame(statistica, risultato), row.names=FALSE) # mostra la tabella dei risultati

#

Dopo avere caricato con la prima riga di codice il pacchetto che include i dati, con la seconda riga di codice al vettore x sono assegnati (<-) i valori della variabile galton$parent che vogliamo utilizzare (si rimanda nuovamente al post Il set di dati galton per la loro storia e il loro significato).

Ora che abbiamo in x i dati che dobbiamo elaborare impieghiamo (terza riga di codice) la funzione length() per calcolarne il numero n.

E finalmente arriviamo (a) al primo loop, che ci serve per calcolare la sommatoria dei dati che chiamiamo sx e che come prima cosa (quarta riga) inizializziamo ponendola uguale a 0 (sx <- 0).

Per calcolare la sommatoria dei dati impieghiamo il loop specificato dalla funzione for() che esegue il codice compreso tra la parentesi graffa aperta { e la parentesi graffa chiusa }.

Il loop for (i in 1:n) procede facendo variare l'indice i da 1 a n e, per ognuno dei valori assunti da i, somma in sx il valore del dato i-esimo x[i]. Alla prima iterazione verrà sommato in sx il valore del dati x₁, alla seconda iterazione verrà aggiunto il valore del dato x₂, e così via fino al dato x_n, dopo il quale il loop si interrompe, fornendo in tal modo la sommatoria dei valori degli n dati. Questa prima parte si conclude con il calcolo della media dei dati ottenuta dividendo la sommatoria dei dati per il loro numero (media <- sx/n).

Il secondo loop (b) calcola la somma dei quadrati delle differenze tra ciascun dato e la media (questa somma viene comunemente denominata devianza). Per questo la variabile devianza che contiene la sommatoria viene prima inizializzata ponendola uguale a 0 (devianza <- 0) quindi con un loop identico al precedente, facendo variare l'indice i da 1 a n, per ognuno dei valori assunti da i viene sommata nella variabile devianza la differenza tra ciascun dato e la media dei dati elevata al quadrato cioè (x[i]-media)^2 ottenendo in tal modo la sommatoria desiderata che divisa per n - 1 ci consente di ottenere la varianza, la cui radice quadrata è infine la deviazione standard, come dalle formule riportate all'inizio.

La terza e ultima parte dello script (c) si limita a generare e a riportare nella Console di R una tabella con la sintesi dei risultati ottenuti:

            statistica  risultato

     numero dei dati = 928.000000

               media =  68.308190

            varianza =   3.194561

 deviazione standard =   1.787333

Il secondo metodo per il calcolo di una sommatoria impiega la vettorializzazione. Molte funzioni in R sono vettorializzate, il che significa che le operazioni specificate all'interno della funzione sono effettuate in parallelo su tutti gli elementi del vettore contenente i dati. 

Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# STATISTICHE ELEMENTARI PARAMETRICHE con vettorializzazione

#

library(psychTools) # carica il pacchetto che include il set di dati galton

x <- galton$parent # salva la variabile nel vettore x

#

n <- length(x) # calcola il numero dei dati

#

### (a)

#

media <- sum(x)/n # calcola la media

#

### (b)

#

varianza <- sum((x-media)^2)/(n-1) # calcola la varianza

ds <- sqrt(varianza) # calcola la deviazione standard

#

### (c)

#

statistica <- c("numero dei dati =", "media =", "varianza =", "deviazione standard =") # predispone i nomi delle statistiche calcolate 

risultato <- c(n, media, varianza, ds) # predispone i risultati delle statistiche calcolate

print(data.frame(statistica, risultato), row.names=FALSE) # mostra la tabella dei risultati

#

Si vede chiaramente in che modo R realizza la vettorializzazione nella quinta riga di codice: la funzione sum() calcola la sommatoria dei quadrati delle differenze tra ciascun dato e la media (la  devianza), necessaria a sua volta per il calcolo della varianza, impiegando direttamente come argomento il vettore x contenente gli n dati, ragion per cui l'espressione (x-media)^2 si intende applicata, e viene applicata, a tutti gli n dati dati del vettore x. Confrontato con quello dello script precedente il codice della parte (a) e della parte (b) è molto più conciso e più facile da leggere, oltre ad essere più efficiente di quello impiegato da un linguaggio non vettorializzato, mentre la parte (c) dello script è semplicemente identica alla precedente (risultati inclusi).

Se volete infine controllare la correttezza dei calcoli effettuati con loop e con vettorializzazione potete confrontarne i risultati con quelli ottenuti con queste tre righe di codice che fanno direttamente ricorso alle funzioni statistiche standard della versione base di R (che si assumono fornire risultati corretti):

#

mean(galton$parent)

var(galton$parent)

sd(galton$parent)

#

In conclusione abbiamo visto attraverso un breve esempio come con R sia possibile effettuare la ripetizione ciclica di un calcolo, o più in generale di una qualsiasi operazione, sia ricorrendo alla classica iterazione (loop), sia impiegando funzioni che prevedono la vettorializzazione [1]. E questo ci ha consentito di illustrare due tecniche di programmazione semplici ma interessanti, che possono risultare utili nel corso dell'utilizzo di R.

----------

[1] R-bloggers. Vectorization in R: Why? Nell'articolo sono riportati anche i riferimenti ad alcune risorse sull'argomento. URL consultato il 20/09/2021: https://bit.ly/3nOsxuS

\bar{x}

Adattare i margini a un grafico

Può capitare di realizzare un grafico per il quale non sono adeguati i margini previsti di default per la finestra grafica di R. Vediamo il problema e una possibile soluzione.

Accertatevi innanzitutto di avere installato il pacchetto DAAG che contiene i dati che ci servono, o in alternativa procedete come indicato nel post Il set di dati ais nel quale trovate anche illustrato il significato dei dati impiegati [1]. Quindi copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# ADATTARE I MARGINI A UN GRAFICO
#
library(DAAG) # carica il pacchetto DAAG incluso il set di dati ais
mydata <- ais[,c(5)] # salva in mydata i valori della variabile contenuta nella colonna 5 della tabella
#
par("mar") # mostra i margini di default della finestra grafica
#
# traccia istogramma e kernel density plot sovrapposti, con i margini di default
#
hist(mydata, main="Istogramma e kernel density plot", xlab="Ferritina in µg/L", ylab = "Numero dei casi per classe", xlim = c(0,250)) # traccia l'istogramma
par(new=TRUE, ann=FALSE) # consente la sovrapposizione del grafico successivo
plot(density(mydata), main="Kernel density plot", xlab="Ferritina in µg/L", ylab = "Frequenza", xlim = c(0,250), yaxt="n", col="red") # traccia il kernel density plot
axis(4, col.ticks="red", col.axis="red") # riporta l'asse delle y sulla destra
mtext("Stima kernel di densità", side=4, line=3) # aggiunge la legenda all'asse delle y sulla destra
#
windows() # apre una nuova finestra grafica
par(mar=c(5.1,4.1,4.1,5.1)) # aumenta il margine destro da 2.1 a 5.1
#
# traccia istogramma e kernel density plot sovrapposti, con margine destro allargato
#
hist(mydata, main="Istogramma e kernel density plot", xlab="Ferritina in µg/L", ylab = "Numero dei casi per classe", xlim = c(0,250)) # traccia l'istogramma
par(new=TRUE, ann=FALSE) # consente la sovrapposizione del grafico successivo
plot(density(mydata), main="Kernel density plot", xlab="Ferritina in µg/L", ylab = "Frequenza", xlim = c(0,250), yaxt="n", col="red") # traccia il kernel density plot
axis(4, col.ticks="red", col.axis="red") # riporta l'asse delle y sulla destra
mtext("Stima kernel di densità", side=4, line=3) # aggiunge la legenda all'asse delle y sulla destra
#
par(mar=c(5.1,4.1,4.1,2.1)) # ripristina i margini di default
#

In questo caso il problema è causato dal desiderio di riportare sullo stesso grafico l'istogramma e il kernel density plot di una variabile. La scala sull'asse delle ascisse x è identica per entrambi, ma nel caso dell'istogramma dobbiamo riportare sull'asse delle ordinate y il numero di casi per ciascuna delle classi in cui è suddiviso, mentre nel caso del kernel density plot dobbiamo riportare sull'asse delle ordinate y la stima kernel di densità espressa in termini di probabilità. La soluzione più ovvia è di riportare due scale delle ordinate, una sulla sinistra riferita all'istogramma, e una sulla destra riferita al kernel density plot. Per chiarezza di rappresentazione riportiamo il kernel density plot in rosso e i relativi valori sulla scale delle ordinate di destra nello stesso colore.

Con la prima parte dello script viene generato questo primo grafico con istogramma e kernel density plot sovrapposti

che impiega i margini di default della finestra grafica, i quali, richiamati nella terza riga di codice mediante la funzione par("mar")

> par("mar") # mostra i margini di default

[1] 5.1 4.1 4.1 2.1

sono rispettivamente 5.1 per il margine inferiore, 4.1 per il margine sinistro, 4.1 per il margine superiore e 2.1 per il margine destro, espressi come numero di righe (per i dettagli digitare help(par) nella Console di R).

Come si vede impiegando la funzione

mtext("Stima kernel di densità", side=4, line=3)

la legenda ("Stima kernel di densità") - aggiunta all'asse delle y sulla destra (side=4) nella terza riga di testo (line=3) - non compare nella finestra grafica.

Nella seconda parte dello script il problema viene corretto impiegando lo stesso identico codice, ma facendolo precedere dalla funzione

par(mar=c(5.1,4.1,4.1,5.1))

che aumenta il margine destro dal valore di default di 2.1 a 5.1 riducendo lievemente la larghezza del grafico e lasciando così sulla destra spazio adeguato per la didascalia dell'asse delle y, che nel grafico precedente non era visibile, e ottenendo questo secondo grafico [2] .

Lo script si conclude con la funzione par(mar=c(5.1,4.1,4.1,2.1)) che ripristina i valori di default dei margini (un passo non obbligatorio, ma altamente raccomandato).

La soluzione riportata può essere facilmente adattata ai diversi casi che si possono presentare nella pratica.

----------

[1] Qui impieghiamo i valori - contenuti nella colonna 5 della tabella ais (vedere la seconda riga di codice dello script) - della variabile ferritina, una proteina che, pur con alcune limitazioni, fornisce una misura dei depositi di ferro presenti nell'organismo, necessari per garantire una adeguata produzione di emoglobina.

[2] Una banalità per chi non fosse ancora abituato a R: quando si esegue lo script, il secondo grafico viene sovrapposto al primo, quindi è necessario minimizzare la finestra del secondo grafico o afferrarla con il mouse e spostarla per vedere il grafico sottostante.

Analisi della varianza a due fattori

Il complesso di tecniche statistiche note come analisi della varianza (ANOVA) è stato introdotto nel post Analisi della varianza a un fattore, nel quale avevamo analizzato i dati di produzione di tre macchine, azionate da operatori non meglio specificati, con la domanda: le differenze tra le medie di produzione delle tre macchine possono essere attribuite al caso? Se così non fosse, dovremmo pensare che esiste qualcosa in grado di influenzare la produzione media delle macchine.

Il disegno sperimentale adottato aveva previsto di raccogliere i dati in modo tale che la variabilità osservata, misurata dalla differenza tra le medie di produzione delle i macchine, risultava essere quella determinata dal fattore macchina+operatore. Questi erano i dati riportati da Wonnacott [1] che abbiamo impiegato

e questo è il grafico che li rappresentava (in nero la media ± 2 deviazioni standard).

Ci proponiamo ora di fare un passo ulteriore: decomporre la variabilità osservata in variabilità dovuta al fattore macchina e variabilità dovuta al fattore operatore. Per questo impieghiamo di nuovo dati illustrati da Wonnacott [2], ma raccolti diversamente, che riportano nelle righe la produzione delle macchine e nelle colonne gli operatori che le azionano. Da notare che le medie [della produzione] delle macchine sono identiche alle precedenti, ma i dati hanno una distribuzione di valori più ampia, un fatto molto importante e che si ripercuote sulle conclusioni che dai dati possiamo trarre.

Il disegno sperimentale adottato ci consente di confrontare tra loro sia le medie delle i macchine (ultima colonna), sia le medie dei j diversi operatori (ultima riga) e analizzare le differenze tra le medie rispondendo separatamente a due domande:

→ le differenze tra le medie delle macchine possono essere attribuite al caso?

→ le differenze tra le medie degli operatori possono essere attribuite al caso? 

Se così non fosse (cioè se le differenze non possono essere attribuite al caso) dovremmo pensare che esiste qualcosa in grado di influenzare la produzione media delle macchine (ad esempio differenze strutturali, inadeguata manutenzione o altro) e/o qualcosa in grado di influenzare la produzione media degli operatori (ad esempio formazione dell'operatore, fenomeni di fatica o altro).

I dati riportati in forma di tabella da Wonnacott per essere analizzati con R devono essere organizzati in un file di testo, sotto forma di righe (record) contenenti ciascuna tre variabili (campi), una prima variabile qualitativa (fattore) che indica la macchina, una seconda variabile qualitativa (fattore) che indica l'operatore e una variabile numerica che indica la produzione della combinazione macchina/operatore, e assumono quindi la forma seguente:

macchina;operatore;produzione
i1;j1;56.7
i1;j2;45.7
i1;j3;48.3
i1;j4;54.6
i1;j5;37.7
i2;j1;64.5
i2;j2;53.4
i2;j3;54.3
i2;j4;57.5
i2;j5;52.3
i3;j1;56.7
i3;j2;50.6
i3;j3;49.5
i3;j4;56.5
i3;j5;44.7

Copiate le sedici righe riportate qui sopra e salvatele in C:\Rdati\ in un file di testo denominato anova2.csv (attenzione all'estensione .csv al momento del salvataggio del file).

In alternativa andate alla pagina Dati nella quale trovate diverse opzioni per scaricare i file di dati, quindi copiate il file anova2.csv nella cartella C:\Rdati\

Poi, se ancora non l'avete fatto, scaricate dal CRAN e installate il pacchetto psych, che ci consentirà di calcolare alcune statistiche riepilogative dei dati, e il pacchetto ggplot2, con il quale realizzeremo una rappresentazione grafica dei dati integrativa all'analisi statistica.

Infine copiate e incollate nella Console di R questo script e premete ↵ Invio:

# ANOVA analisi della varianza a due fattori

#

mydata <- read.table("c:/Rdati/anova2.csv", header=TRUE, sep=";", dec=".") # importa i dati [1]

anova2 <- aov(produzione~macchina+operatore, data=mydata) # calcola la variabilità tra macchine e tra operatori

summary(anova2) # mostra i risultati

#

library(psych) # carica il pacchetto

describeBy(mydata$produzione, mydata$macchina) # statistiche della produzione per macchina 

describeBy(mydata$produzione, mydata$operatore) # statistiche della produzione per operatore

#

pairwise.t.test(mydata$produzione, mydata$macchina, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie delle macchine impiegando la correzione di bonferroni per confronti multipli

pairwise.t.test(mydata$produzione, mydata$operatore, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie degli operatori impiegando la correzione di bonferroni per confronti multipli [2]

#

library(ggplot2) # carica il pacchetto

plot1 <- ggplot(mydata, aes(x=macchina, y=produzione, fill=macchina)) + geom_dotplot(binaxis='y', stackdir='center', stackratio=1, dotsize=1, binwidth=0.4, show.legend=FALSE) + coord_cartesian(ylim=c(40, 70)) + labs(title="Variabilità tra le macchine", x="Macchina impiegata", y="Produzione realizzata") + theme_classic() # dotplot della produzione per macchina

plot2 <- ggplot(mydata, aes(x=operatore, y=produzione, fill=operatore)) + geom_dotplot(binaxis='y', stackdir='center', stackratio=1, dotsize=1, binwidth=0.4, show.legend=FALSE) + coord_cartesian(ylim=c(40, 70)) + labs(title="Variabilità tra gli operatori", x="Operatore", y="Produzione realizzata") + theme_classic() # dotplot della produzione per operatore

library(gridExtra) # carica il pacchetto

grid.arrange(plot1, plot2, nrow = 1) # i due dotplot sono mostrati affiancati orizzontalmente

#

Lo script è suddiviso in quattro blocchi di codice. Nel primo blocco, dopo avere importato i dati con la funzione read.table() (prima riga), viene eseguita l'analisi della varianza a un fattore mediante la funzione aov() (seconda riga) e sono riepilogati (terza riga) i risultati con la funzione summary():

            Df Sum Sq Mean Sq F value   Pr(>F)    

macchina     2  154.8   77.40   13.10 0.002997 ** 

operatore    4  381.7   95.43   16.15 0.000673 ***

Residuals    8   47.3    5.91                     

---

Signif. Codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

La varianza (Mean Sq) è calcolata dividendo la somma degli scarti quadratici (Sum Sq) per i gradi di libertà (Df). Il valore F (F value) viene calcolato separatamente per la macchina e per l'operatore, dividendo la varianza  dovuta alle macchine (77.40) per la varianza residua (5.91) ovvero dividendo la varianza  dovuta agli operatori (95.43) per la varianza residua (5.91) .

La probabilità (Pr(>F)) è la probabilità di osservare per caso il valore F e consente di rispondere alla domanda: le differenze tra le medie di produzione possono essere attribuite al caso? Per le medie di produzione delle macchine abbiamo un valore Pr(>F)= 0.002997  e per le medie di produzione degli operatori abbiamo un valore Pr(>F)= 0.000673. Poiché entrambi sono inferiori al valore 0.05 assunto in genere come valore soglia, concludiamo che le differenze osservate non sono attribuibili al caso o, se preferite, che sono significative: quindi esiste qualcosa in grado di influenzare la produzione a livello delle macchine e/o a livello degli operatori.

Il secondo blocco di codice prevede di caricare il pacchetto psych e impiegare la funzione describeBy() per riportare le statistiche elementari prima per le tre macchine i e poi per i cinque operatori j. Si tratta di un passo importante, che ci mostra che i dati, come richiesto dagli assunti alla base dell'ANOVA [3], sono distribuiti in modo gaussiano (al bisogno per l'interpretazione si rimanda all'analisi elementare di dati univariati riportata nel post Analisi esplorativa dei dati).

> describeBy(mydata$produzione, mydata$macchina) # statistiche della produzione per macchina 

 Descriptive statistics by group 

group: i1

   vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se

X1    1 5 48.6 7.57   48.3    48.6 9.34 37.7 56.7    19 -0.26    -1.79 3.38

------------------------------------------------------------ 

group: i2

   vars n mean   sd median trimmed  mad  min  max range skew kurtosis  se

X1    1 5 56.4 4.93   54.3    56.4 2.97 52.3 64.5  12.2 0.72    -1.41 2.2

------------------------------------------------------------ 

group: i3

   vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se

X1    1 5 51.6 5.08   50.6    51.6 8.75 44.7 56.7    12 -0.14    -1.93 2.27

> describeBy(mydata$produzione, mydata$operatore) # statistiche della produzione per operatore

 Descriptive statistics by group 

group: j1

   vars n mean  sd median trimmed mad  min  max range skew kurtosis  se

X1    1 3 59.3 4.5   56.7    59.3   0 56.7 64.5   7.8 0.38    -2.33 2.6

------------------------------------------------------------ 

group: j2

   vars n mean  sd median trimmed  mad  min  max range  skew kurtosis   se

X1    1 3 49.9 3.9   50.6    49.9 4.15 45.7 53.4   7.7 -0.17    -2.33 2.25

------------------------------------------------------------ 

group: j3

   vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se

X1    1 3 50.7 3.17   49.5    50.7 1.78 48.3 54.3     6 0.32    -2.33 1.83

------------------------------------------------------------ 

group: j4

   vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se

X1    1 3 56.2 1.47   56.5    56.2 1.48 54.6 57.5   2.9 -0.2    -2.33 0.85

------------------------------------------------------------ 

group: j5

   vars n mean  sd median trimmed   mad  min  max range skew kurtosis   se

X1    1 3 44.9 7.3   44.7    44.9 10.38 37.7 52.3  14.6 0.03    -2.33 4.22

Il terzo blocco di codice consiste in due righe di codice, una per analizzare i dati delle macchine, e una analizzare i dati degli operatori. Abbiamo già visto che l'analisi della varianza a un fattore è una generalizzazione del test t di Student per dati non appaiati, ed è equivalente al test t quando i gruppi sono due (vedere il post Test parametrici e non parametrici per due campioni indipendenti). Tuttavia esistono delle correzioni del test t che consentono di impiegarlo anche nel caso del confronto di più di due gruppi (vedere il post Test parametrici e non parametrici per più campioni indipendenti). Queste correzioni controbilanciano l'aumento della probabilità, derivante dall'esecuzione di confronti multipli, di considerare la differenza tra le medie di due campioni come significativa quando invece non è significativa, ed è proprio una di queste, la correzione di Bonferroni (p.adjust.method = "bonferroni"), che viene qui applicata ai dati mediante la funzione pairwise.t.test() [3]. Mentre l'ANOVA è un test globale, che ci dice che tra le medie esiste una qualche differenza, ma non ci consente di individuare la/e media/e responsabile/i di tale differenza, il test t viene effettuato per tutti i confronti possibili tra medie.

Il test t effettuato per tutti i confronti possibili tra macchine indica differenze sempre non significative - essendo sempre superiori a 0.05 le probabilità che le differenze tra medie qui osservate siano dovute al caso:

> pairwise.t.test(mydata$produzione, mydata$macchina, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie delle macchine impiegando la correzione di bonferroni per confronti multipli

        Pairwise comparisons using t tests with pooled SD 

data:  mydata$produzione and mydata$macchina 

   i1   i2  

i2 0.18 -   

i3 1.00 0.69

P value adjustment method: bonferroni 

Il test t effettuato per tutti i confronti possibili tra operatori ci indica che la significatività dall'ANOVA è dovuta ad un'unica differenza tra tutti i confronti effettuati: j1 risulta diverso da j5 con un valore p=0.029, di poco inferiore al valore soglia 0.05 comunemente adottato:

> pairwise.t.test(mydata$produzione, mydata$operatore, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie degli operatori impiegando la correzione di bonferroni per confronti multipli

        Pairwise comparisons using t tests with pooled SD 

data:  mydata$produzione and mydata$operatore 

   j1    j2    j3    j4   

j2 0.283 -     -     -    

j3 0.411 1.000 -     -    

j4 1.000 1.000 1.000 -    

j5 0.029 1.000 1.000 0.117

P value adjustment method: bonferroni 

Quindi il test t con la correzione di Bonferroni:

→ nel confronto tra le medie delle macchine indica differenze sempre non significative, contrariamente all'ANOVA;

→ nel confronto tra le medie degli operatori conferma il risultato dell'ANOVA specificando che la differenza è imputabile ad un unico caso, quello degli operatori j1 e j5.

Il quarto e ultimo blocco di codice prevede di caricare la libreria ggplot2 per realizzare un grafico a punti (dotplot) che visualizza i dati della variabilità tra le macchine e un grafico a punti che visualizza i dati delle variabilità tra gli operatori, e combinarli in un'unica figura.

Per i dettagli delle funzioni e degli argomenti impiegati per questa rappresentazione grafica si rimanda al post Grafici a punti (dotplot) mentre nel post Grafici a violino (violin plot) trovate anche come sovrapporre ai punti un grafico a violino o un grafico a scatola con i baffi (boxplot).

Se l'analisi della varianza consente di effettuare un confronto tra medie, a questo punto sembra logico sovrapporre ai dati raccolti la loro media e la loro deviazione standard.  Lo si può fare aggiungendo di seguito all'ultima riga di codice dello script riportato sopra subito dopo + theme_classic() questo codice

+ stat_summary(fun.data=mean_sdl, fun.args = list(mult=2), geom="pointrange", color="black", show.legend=FALSE) # aggiunge media +/- 2 deviazioni standard

Da notare che nella funzione stat_summary() l'argomento

mult=2

specifica di rappresentare attorno alla media l'intervallo uguale a 2 deviazioni standard, ma potete anche mettere

mult=1

per mostrare l'intervallo uguale a una deviazione standard, o qualsiasi altro valore per l'intervallo che desiderate rappresentare.

In alternativa copiate e incollate nella Console di R queste due righe di codice che aprono una nuova finestra grafica nella quale viene rappresentato il grafico che sovrappone ai punti la media con l'intervallo corrispondente a due deviazioni standard:

#
windows() # apre una nuova finestra grafica
library(ggplot2) # carica il pacchetto
plot1 <- ggplot(mydata, aes(x=macchina, y=produzione, fill=macchina)) + geom_dotplot(binaxis='y', stackdir='center', stackratio=1, dotsize=1, binwidth=0.4, show.legend=FALSE) + coord_cartesian(ylim=c(40, 70)) + labs(title="Variabilità tra le macchine", x="Macchina impiegata", y="Produzione realizzata") + theme_classic() + stat_summary(fun.data=mean_sdl, fun.args = list(mult=2), geom="pointrange", color="black", show.legend=FALSE) # dotplot della produzione per macchina
plot2 <- ggplot(mydata, aes(x=operatore, y=produzione, fill=operatore)) + geom_dotplot(binaxis='y', stackdir='center', stackratio=1, dotsize=1, binwidth=0.4, show.legend=FALSE) + coord_cartesian(ylim=c(40, 70)) + labs(title="Variabilità tra gli operatori", x="Operatore", y="Produzione realizzata") + theme_classic() + stat_summary(fun.data=mean_sdl, fun.args = list(mult=2), geom="pointrange", color="black", show.legend=FALSE) # dotplot della produzione per macchina
#

library(gridExtra) # carica il pacchetto

grid.arrange(plot1, plot2, nrow = 1) # i due dotplot sono mostrati affiancati orizzontalmente

#

La rappresentazione grafica dei dati ci viene nuovamente in aiuto nell'interpretazione dei risultati. Come si vede in questo caso i dati sono molto dispersi, e gli intervalli media ± 2 deviazioni standard sono ampiamente sovrapposti: e questo corrobora la non significatività delle differenze tra le medie. Da notare la notevole differenza rispetto ai dati dell'ANOVA a un fattore riportati nella prima figura in alto. La maggior differenza rilevata graficamente è quella tra operatore j1 e operatore j5 ed è in linea con la debole (appena inferiore a 0.05) significatività statistica (p=0.029) di tale differenza. 

La conclusione? Impiegando i dati tratti da un testo di statistica abbiamo sviluppato due esempi di analisi della varianza, a un fattore [post Analisi della varianza a un fattore] e a due fattori [questo post], dai quali possiamo ricavare alcune considerazioni interessanti:

→ l'analisi della varianza (ANOVA) a dispetto del nome, è un insieme di tecniche per effettuare confronti multipli tra medie;

→ l'ANOVA fornisce come risultato un rapporto tra varianze (il test F), che è un test globale che ci dice che tra le medie esiste una qualche differenza, ma non specifica la/e media/e responsabile/i di tale differenza;

→ nel caso limite in cui il confronto tra medie è limitato a due campioni l'ANOVA equivale al test t di Student per dati non appaiati (campioni indipendenti);

→ i confronti multipli tra medie possono essere effettuati anche impiegando test alternativi come il test t con opportune correzioni, ad esempio con la correzione di Bonferroni, con il vantaggio, in questo caso, che la significatività viene valutata separatamente per ciascuna coppia di medie poste a confronto;

→ il valore soglia p=0.05 non è un dogma, e qualora si ritenga opportuno essere più prudenti (conservativi) nel giudicare la significatività di un test statistico si può adottare un valore soglia inferiore come per esempio p=0.01;

→ quando i dati sono poco dispersi e l'intervallo media ± 2 deviazioni standard non si sovrappone (vedere la prima immagine riportata qui sopra e il post Analisi della varianza a un fattore) i risultati di ANOVA e test t (con le opportune correzioni) sono uguali in quanto l'informazione fornita dai dati è elevata, lascia poco adito a dubbi [sulla significatività delle differenze tra le medie] e la diversità degli assunti alla base delle tecniche statistiche impiegate non porta a conclusioni diverse;

→ quando i dati sono molto dispersi e l'intervallo media ± 2 deviazioni standard è ampiamente sovrapposto (vedere la seconda e la terza immagine riportate qui sopra), i risultati di ANOVA e test t (con le opportune correzioni) possono differire in quanto l'informazione fornita dai dati è scarsa, lascia adito a molti dubbi [sulla significatività delle differenze tra le medie] e la diversità degli assunti alla base delle tecniche statistiche impiegate può portare a conclusioni diverse;

→ la rappresentazione grafica dei dati rappresenta come sempre una importante integrazione ai test statistici.

In sintesi anche i risultati dell'ANOVA, come del resto tutti i risultati dell'analisi statistica, non devono essere interpretati in modo schematico e rigido, ma devono essere interpretati:

→ verificando che i dati analizzati rispettino gli assunti dell'ANOVA;

→ integrando i risultati dell'ANOVA con quelli di test alternativi per il confronto tra medie;

→ valutando criticamente le eventuali discrepanze tra i risultati dell'ANOVA e dei test alternativi;

→ adottando al bisogno requisiti di significatività più stringenti del tradizionale p=0.05;

→ integrando i risultati dell'ANOVA e dei test alternativi con l'esplorazione grafica dei dati;

→ traendo le conclusioni sulla base di una valutazione globale dei risultati ottenuti.

Il che ci ricorda che tutto sommato la statistica non dovrebbe essere altro che buonsenso, reso scientifico nella misura in cui si trasforma in un modello quantitativo e numerico, e adiuvato dalla rappresentazione grafica.

----------

[1] Wonnacott TH, Wonnacott RJ. Introduzione alla statistica. Franco Angeli Editore, Milano, 1980, ISBN 88-204-0323-4, Tabella 10-1, p. 238.

[2] Wonnacott TH, Wonnacott RJ. Introduzione alla statistica. Franco Angeli Editore, Milano, 1980, ISBN 88-204-0323-4, Tabella 10-9, p. 255.

[3] L'analisi della varianza è basata su due ipotesi: (i) che i dati siano distribuiti in modo gaussiano e (ii) che la varianza sia la stessa nei diversi gruppi confrontati.