Analisi della varianza a due fattori

Il complesso di tecniche statistiche note come analisi della varianza (ANOVA) è stato introdotto con l'analisi della varianza a un fattore con la 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) e che mostra come medie risultassero tutte e tre significativamente diverse l'una dall'altra.

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 aggiungendo un ↵ Invio al termine dell'ultima riga 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\

Scaricate dal CRAN e installate 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

#

pairwise.t.test(mydata$produzione, mydata$macchina, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie delle macchine impiegando il test t con 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 il test t con la correzione di bonferroni per confronti multipli [2]

#

Lo script è suddiviso in tre blocchi di codice. Nel primo blocco sono importati i dati con la funzione read.table().

Nel secondo 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 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 [4]. Tuttavia esistono delle correzioni del test t che consentono di impiegarlo anche nel caso del confronto di più di due gruppi [5]. 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(). 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 mentre l'ANOVA indica la presenza nei dati di una qualche differenza significativa, il test t con la correzione di Bonferroni consente di rilevare dove questa si trova in quanto:

→ nel confronto tra le medie delle macchine indica differenze sempre non significative;

→ nel confronto tra le medie degli operatori indica che la differenza è imputabile ad un unico caso, quello degli operatori j1 e j5.

Se oltre alla libreria ggplot2 installate la libreria gridExtra potete realizzare un grafico a punti [6] 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.

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

#

library(ggplot2) # carica il pacchetto per la grafica

library(gridExtra) # carica il pacchetto per combinare i grafici

#

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

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

#

Se l'ANOVA consente di effettuare un confronto tra medie, a questo punto sembra logico sovrapporre ai dati raccolti la loro media e la loro deviazione standard.

Copiate e incollate nella Console di R queste quattro 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 e inizializza una nuova finestra grafica

#

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 +/- 2 deviazioni standard

#

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 +/- 2 deviazioni standard

#

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

#

Da notare che nella funzione stat_summary() l'argomento mult=2 che specifica di rappresentare attorno alla media l'intervallo uguale a 2 deviazioni standard può essere modificato.

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 [7] e a due fattori, 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;

→ perché l'ANOVA fornisca risultati affidabili è necessario verificare in via preliminare che i dati siano distribuiti in modo gaussiano e con varianze omogenee, eseguendo i test opportuni;

→ 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 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 di normalità e omogeneità delle varianze previsti dall'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 in statistica il rigore scientifico dei modelli matematici e dei numeri dovrebbe sempre essere integrato con il buonsenso, al quale può contribuire in modo importante una adeguata 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.

[4] Vedere il post Test parametrici e non parametrici per due campioni indipendenti.

[5] Vedere il post Test parametrici e non parametrici per più campioni indipendenti.

[6] 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).

Analisi della varianza a un fattore

Anche se può sembrare paradossale, l'analisi della varianza consente di effettuare il confronto tra medie. Infatti come è stato sottolineato, a dispetto della “... apparente incoerenza della nomenclatura – il fatto cioè che le tecniche riguardino confronti di medie piuttosto che varianze – ... il complesso delle tecniche chiamato analisi della varianza costituisce un potente metodo per analizzare il modo in cui il valore medio di una variabile è influenzato da classificazioni di vario genere dei dati” [1].

Ancor più direttamente Wonnacott all'inizio del capitolo dedicato all'analisi della varianza riporta: "... abbiamo compiuto inferenze relative alla media di una popolazione ... abbiamo esteso il ragionamento alla differenza fra due medie ... vogliamo confrontare ora r medie, usando un insieme di tecniche che vengono comunemente denominate «analisi della varianza»" [2].  

Le tecniche di analisi della varianza (ANOVA) sono ampiamente trattate in tutti i testi di statistica, ai quali si rimanda per i necessari approfondimenti. Vediamo in questo post la più semplice di queste tecniche, l'analisi della varianza a un fattore che:

→ è "... una generalizzazione del test t [di Student] per [il confronto tra medie di] dati non appaiati [due campioni indipendenti], adatta a un numero qualunque di gruppi ..." [1];

→ è "... equivalente al test t per [il confronto tra medie nel caso di] dati non appaiati quando i gruppi sono ... due" [1];

→ è basata "... su due importanti ipotesi: (a) la normalità delle distribuzioni delle osservazioni ... e (b) la costanza delle varianze nei diversi gruppi" [3].

Per la verifica preliminare della omogeneità delle varianze (b) si rimanda al post che tratta il test di Levene nel quale sono riportati anche i riferimenti ai test che possono essere eseguiti per verificare (a) ovvero la normalità delle distribuzioni dei dati [4].

I problemi che si presentano nella pratica quando gli assunti alla base del modello statistico dell'ANOVA non sono completamente soddisfatti – e i correttivi che possono eventualmente essere adottati – vanno al di la dei limiti di questo blog, ma sono ben approfonditi per esempio da Snedecor [5].

Per lo script impieghiamo i dati riportati da Wonnacott che riguardano "... tre macchine (A, B e C), le quali, essendo azionate da uomini e a causa di altre ragioni inesplicabili, danno luogo ad un prodotto orario soggetto a fluttuazioni casuali. Nella speranza di «mediare» e quindi di ridurre gli effetti di tali fluttuazioni, si effettua un campione casuale di 5 ore per ciascuna macchina, i cui risultati sono raccolti nella Tabella 10-1 insieme con le relative medie" [6].

La domanda è: le differenze tra le medie di produzione 48.6, 56.4 e 51.6, riportate nell'ultima colonna sulla destra, possono essere attribuite al caso? Se così non fosse, dovremmo pensare che esiste qualcosa in grado di influenzare la produzione media delle macchine.

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 due variabili (campi), una variabile qualitativa (fattore) che indica la macchina, e una variabile numerica che indica la produzione della macchina, e assumono quindi la forma seguente:

macchina;produzione
i1;48.4
i1;49.7
i1;48.7
i1;48.5
i1;47.7
i2;56.1
i2;56.3
i2;56.9
i2;57.6
i2;55.1
i3;52.1
i3;51.1
i3;51.6
i3;52.1
i3;51.1

Copiate le sedici righe riportate qui sopra aggiungendo un ↵ Invio al termine dell'ultima riga e salvatele in C:\Rdati\ in un file di testo denominato anova1.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 anova1.csv nella cartella C:\Rdati\

Poi scaricate dal CRAN e installate 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 un fattore
#

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

#

anova1 <- aov(produzione~macchina, data=mydata) # esegue l'analisi della varianza
summary(anova1) # mostra i risultati dell'analisi della varianza
#
pairwise.t.test(mydata$produzione, mydata$macchina, p.adjust.method = "bonferroni") # confronto tra le (produzioni) medie delle macchine impiegando il test t con la correzione di Bonferroni per confronti multipli
#

Lo script è suddiviso in tre blocchi di codice. Nel primo sono importati i dati con la funzione read.table().

Nel secondo 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():

> summary(anova1) # mostra i risultati dell'analisi della varianza

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

macchina     2 154.80   77.40   141.6 4.51e-09 ***
Residuals   12   6.56    0.55                     
---
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 (141.6) del rapporto fra varianze (F value) viene calcolato dividendo la varianza dovuta alle macchine (77.40) per la varianza residua (0.55).

La probabilità (Pr(>F)) è la probabilità di osservare per caso il valore F (141.6) e consente di rispondere alla domanda: le differenze tra le medie di produzione possono essere attribuite al caso? Essendo tale probabilità 4.51e-09 o se preferite 0.0000000451, quindi molto inferiore al valore 0.05 assunto in genere come valore soglia, concludiamo che le differenze tra le medie di produzione non sono attribuibili al caso: le medie sono significativamente diverse, quindi esiste qualcosa in grado di influenzare la produzione [media] delle macchine. 

Da notare una cosa molto importante: l'ANOVA è un test globale e ci dice che tra le medie [nel nostro caso tra le medie di produzione delle tre macchine] esiste una differenza significativa, ma non consente di individuare la/e media/e causa della significatività. In altre parole non ci dice se la significatività sia dovuta alla differenza tra la prima macchina e la seconda, tra la prima macchina e la terza, tra la seconda macchina e la terza, o a una qualche combinazione di queste tre possibilità. Se è questo che interessa, è possibile ricorrere a test alternativi (vedi sotto).

Il terzo blocco di codice consiste anch'esso in una sola riga. 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 [7]. Tuttavia esistono delle correzioni del test t che consentono di impiegarlo anche nel caso del confronto di più di due gruppi [8], 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() per confrontarne i risultati con quelli dell'ANOVA

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 3.4e-09 -      

i3 1.0e-04 8.1e-07

P value adjustment method: bonferroni 

Come si vede il test t effettuato per tutti e tre i confronti possibili tra macchine conferma i risultati dell'ANOVA riportando differenze sempre significative – essendo sempre inferiori a 0.05 le probabilità che le differenze tra medie qui osservate siano dovute al caso – ma ci indica anche dove risiedono tali differenze, riportando separatamente che i1 <> i2 (p=3.4e-09), i1 <> i3 (p=1.0e-04), i2 <> i3 (p=8.1e-07), cosa che l'ANOVA, essendo un test globale, non consente di specificare (vedi sopra). 

Altre correzioni alternative, meno conservative della correzione di Bonferroni [9], possono essere impiegate con l'argomento p.adjust.method = : secondo Holm (1979) ("holm"), Hochberg (1988) ("hochberg"), Hommel (1988) ("hommel"), Benjamini & Hochberg (1995) ("BH" o "fdr"), Benjamini & Yekutieli (2001) ("BY"), e "none" per nessuna correzione.

Se installate  la libreria ggplot2 potete realizzare un grafico a punti (dotplot) [10] che consente di visualizzare i dati: il che è sempre un importante complemento all'analisi numerica. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

 

#

library(ggplot2) # carica il pacchetto per la grafica

#

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, 60)) + labs(title="Variabilità osservata nella produzione tra le macchine", x="Macchina impiegata", y="Produzione realizzata") + theme_classic() # dotplot della produzione per macchina

#

Se l'ANOVA consente di effettuare un confronto tra medie, a questo punto sembra logico sovrapporre ai dati raccolti la loro media e la loro deviazione standard.

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

#

windows() # apre e inizializza una nuova finestra grafica

#
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, 60)) + labs(title="Variabilità osservata nella produzione 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) # aggiunge media +/- 2 deviazioni standard

#

Da notare che nella funzione stat_summary l'argomento mult=2 che specifica di rappresentare attorno alla media l'intervallo uguale a 2 deviazioni standard può essere modificato. 

La rappresentazione grafica dei dati ci viene nuovamente in aiuto nell'interpretazione dei risultati: come si vede i dati sono poco dispersi, le medie sono ben differenti, gli intervalli media ± 2 deviazioni standard non si sovrappongono, e questo conferma la significatività delle differenze tra le medie.

Dal punto di vista statistico è interessante notare l'importanza del disegno sperimentale. Quello qui adottato prevedeva di raccogliere i dati in modo tale che la variabilità osservata è quella dovuta al fattore macchina+operatore, pertanto non è possibile stabilire se le differenze (tra medie) osservate sono imputabili alla sola macchina, al solo operatore o a entrambi.

Impiegando un altro disegno sperimentale, e raccogliendo i dati diversamente, è possibile decomporre la variabilità osservata in variabilità dovuta al fattore macchina e variabilità dovuta al fattore operatore. lo facciamo nel post Analisi della varianza a due fattori, al termine del quale trovate anche alcune considerazioni generali sull'ANOVA.

----------

[1] Armitage P. Statistica medica. Giangiacomo Feltrinelli Editore, Milano, 1979, p. 188.

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

[3] Armitage P. Statistica medica. Giangiacomo Feltrinelli Editore, Milano, 1979, pp. 195-196.

[4] Vedere il post Come valutare l'omogeneità tra varianze.

[5] Snedecor GW, Cochran WG. Statistical Methods. The Iowa State University Press, 1980, ISBN 0-8138-1560-6. Chapter 15 - Failures in the assumptions, pp. 274-297.

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

[7] Vedere il post Test parametrici e non parametrici per due campioni indipendenti.

[8] Vedere il post Test parametrici e non parametrici per più campioni indipendenti per il significato di queste correzioni, che 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.

[9] Un test statistico più conservativo a parità di condizioni produce come risultato un valore di p più alto quindi ci aspettiamo che a lungo andare fornisca un numero inferiore di differenze significative, mentre un test statistico meno conservativo a parità di condizioni produce come risultato un valore di p più basso quindi ci aspettiamo che a lungo andare fornisca un numero maggiore di differenze significative.

[10] Si rimanda al post Grafici a punti (dotplot), inoltre 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).

Regressione lineare semplice parametrica

La regressione lineare tout court, quella illustrata in tutti i testi di statistica e qui sviluppata con R, è:

→ una regressione lineare semplice (in contrapposizione alla regressione lineare multipla);

→ parametrica (in contrapposizione alla regressione lineare non parametrica);

→ x variabile indipendente (in contrapposizione ad alternative che non prevedono questo assunti).

Dato che la denominazione razionale, quella completa, diverrebbe chilometrica, e dato che l'espressione regressione lineare è però troppo generica, laddove è opportuno semplificare impiegherò la denominazione regressione lineare ordinaria, con riferimento al fatto che è quella non solo dovunque illustrata ma anche più ampiamente impiegata.

Gli assunti alla base del modello matematico-statistico implicano una serie di requisiti che devono essere soddisfatti dai dati, che sono ben chiariti ad esempio in Marubini [1] e in Snedecor [2], ma anche online [3]. Per il caso speciale, ma non infrequente, nel quale nessuna delle due variabili confrontate abbia i requisiti richiesti per essere considerata come variabile indipendente, vedere anche [4].

Qualora invece sia necessario l'impiego di "metodi robusti" di regressione - cioè di metodi che risentono poco della eventuale presenza di dati apparentemente anomali ma che non possono essere esclusi a priori dal calcolo della regressione - è possibile impiegare un metodo non parametrico [5].

Tornando alla nostra regressione lineare semplice parametrica x variabile indipendente, essendo x la variabile indipendente posta sull'asse delle ascisse, y la variabile dipendente posta sull'asse delle ordinate, ed essendo soddisfatti i requisiti previsti per i dati, il metodo dei minimi quadrati consente di calcolare l'intercetta a e il coefficiente angolare b della retta di regressione di equazione

y = a + b∙x

che meglio approssima la distribuzione dei dati sperimentali.

Il calcolo dell'equazione della retta di regressione viene effettuato mediante la funzione lm(), che può essere applicata anche al caso di più variabili indipendenti, consentendo quindi il calcolo della regressione lineare multipla [4]. Ma la cosa più interessante è che R, oltre al calcolo dell'equazione della retta di regressione, un calcolo di per sé semplice, fornisce una serie molto interessante di strumenti che consentono di valutare quanto i dati soddisfano i requisiti richiesti, o in altre parole di valutare se la regressione lineare descrive in modo adeguato la relazione tra le due variabili.

Qui impieghiamo come esempi due diversi set di dati, per uno solo dei quali, come vedremo, la regressione lineare fornisce risultati adeguati.

Il primo è il set di dati ais, nel quale prendiamo in considerazione la concentrazione degli eritrociti (espressa in 10¹²/L) e il valore ematòcrito (espresso in %), due variabili che, analizzate mediante i coefficienti di correlazione parametrici e non parametrici, hanno mostrato di essere correlate, e che, all'ispezione visiva di una serie di grafici di dispersione (xy) realizzati mediante i correlogrammi, hanno mostrato valori ben allineati su una possibile retta.

Per procedere dovete scaricare e installare dal CRAN i pacchetti aggiuntivi gvlma e car, oltre al pacchetto aggiuntivo DAAG [6].

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

# REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x

#

library(DAAG) # carica il pacchetto che include il set di dati ais

library(gvlma) # carica il pacchetto per la funzione gvlma()

library(car) # carica il pacchetto per analisi ouliers e analisi grafica

str(ais) # mostra la struttura dei dati

#

var_x <- ais$rcc # variabile eritrociti x in ascisse

var_y <- ais$hc # variabile ematòcrito y in ordinate

reglin <- lm(var_y ~ var_x) # calcola intercetta (a) e coefficiente angolare (b)

coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc

confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare

#

# analisi statistica della adeguatezza della regressione lineare

#

summary(reglin) # mostra un riepilogo dei risultati

t.test(residuals(reglin)) # verifica che la media degli errori non sia significativamente diversa da zero

shapiro.test(residuals(reglin)) # verifica la normalità della distribuzione degli errori

summary(gvlma(reglin)) # test globale per l'assunto di linearità

outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)

#

# analisi grafica della adeguatezza della regressione lineare

#

windows() # apre una nuova finestra

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

newx = seq(min(var_x), max(var_x), by = 0.01) # valori della x per i quali calcolare l'intervallo di confidenza

conf_interval <- predict(reglin, newdata=data.frame(var_x=newx), interval="confidence", level = 0.95) # calcola gli intervalli di confidenza

plot(var_x, var_y, xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Regressione lineare y = a + b·x") # grafico dei dati

abline(reglin, col="lightblue") # retta di regressione

lines(newx, conf_interval[,2], col="blue", lty=2) # limite di confidenza inferiore

lines(newx, conf_interval[,3], col="blue", lty=2) # limite di confidenza superiore

#

plot(var_y, residuals(reglin), xlab="Ematòcrito (%)", ylab="Ematòcrito osservato - calcolato (%)", main="Analisi delle differenza residue") # grafico delle differenza tra ematòcrito osservato e ematòcrito calcolato con l'equazione della retta

#

influencePlot(reglin, fill=FALSE, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

#

qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

#

Dopo avere caricato i pacchetti aggiuntivi e i dati, e dopo avere mostrato con str(ais) la struttura di questi ultimi, con var_x <- ais$rcc e con var_y <- ais$hc sono memorizzate negli oggetti var_x e var_y rispettivamente la variabile indipendente x, posta in ascisse, e la variabile dipendente y, posta in ordinate. In questo modo sarà possibile riutilizzare per intero lo script semplicemente sostituendo ad ais$rcct e ad ais$hc i vettori contenenti i propri dati.

Con reglin <- lm(var_y ~ var_x) l'equazione della retta di regressione che esprime la y in funzione della x viene calcolata e memorizzata nell'oggetto reglin, che a questo punto diventa l'argomento chiave, l'argomento contenente i dati in ingresso impiegati nelle funzioni successive [7].

Con le funzioni coefficients() e confint() sono infine mostrati i coefficienti della retta di regressione e i loro intervalli di confidenza al 95%

> coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc 

(Intercept)       var_x 

   8.183033    7.398052 

> confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare 

               2.5 %    97.5 %

(Intercept) 6.173717 10.192350

var_x       6.974206  7.821898

e pertanto questa risulta essere l'equazione della retta di regressione:

hc = 8.183033 + 7.398052 · rcc

A questo punto seguono due blocchi di codice, il primo per effettuare una analisi statistica, e il secondo per effettuare una analisi grafica della regressione, entrambe allo scopo di valutare, come già detto, se la regressione lineare descrive in modo adeguato la relazione tra le due variabili.

Per quanto concerne l'analisi statistica, le principali conclusioni sono quelle tratte con il test globale per l'assunto di linearità effettuato mediante la funzione gvlma(), che conferma il fatto che gli assunti che stanno alla base del modello di regressione lineare sono tutti soddisfatti:

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                      Value p-value                Decision

Global Stat        3.479904  0.4809 Assumptions acceptable.

Skewness           2.233166  0.1351 Assumptions acceptable.

Kurtosis           0.958154  0.3277 Assumptions acceptable.

Link Function      0.005239  0.9423 Assumptions acceptable.

Heteroscedasticity 0.283345  0.5945 Assumptions acceptable.

Il test di Bonferroni non evidenzia dati anomali, ma segnala il dato numero 68 come quello che maggiormente si discosta dai rimanenti:

> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers) 

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

    rstudent unadjusted p-value Bonferonni p

68 -3.569525         0.00044828     0.090552

Le conclusioni dell'analisi grafica confermano, con i valori delle differenze residue dispersi in modo casuale, che la regressione lineare descrive in modo adeguato la relazione tra le due variabili

e confermano anche la presenza di possibili dati anomali con il grafico delle distanze di Cook e con il grafico dei quantili, dati che sono evidenziati sia direttamente nei grafici, sia nei rispettivi riepiloghi riportati nella Console di R:

> influencePlot(reglin, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni 

      StudRes         Hat      CookD

68  -3.569525 0.006301040 0.03815682

78   3.186317 0.009724284 0.04766681

161 -2.179789 0.039758811 0.09655642

166  1.363972 0.099962743 0.10287127

> #  

> qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati 

[1] 68 78

I dati per i quali i metodi di identificazione meglio concordano sono il numero 68 e il numero 78, ma anche i dati numero 161 e numero 166 meritano di essere valutati. Si tratta di casi dei quali sarebbe importante controllare la validità, o che potrebbero essere situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati.

Vediamo ora la stessa identica analisi applicata al set di dati galton [8]. I pacchetti aggiuntivi psychTools, gvlma e car, se non l'avete già fatto, devono essere preventivamente scaricati e installati dal CRAN [9].

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

# REGRESSIONE LINEARE SEMPLICE PARAMETRICA y = a + b · x

#

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

library(gvlma) # carica il pacchetto per la funzione gvlma()

library(car) # carica il pacchetto per analisi outliers e analisi grafica

str(galton) # mostra la struttura dei dati

#

var_x <- galton$parent # variabile altezza dei padri x in ascisse

var_y <- galton$child # variabile altezza dei figli y in ordinate

reglin <- lm(var_y ~ var_x) # calcola intercetta (a) e coefficiente angolare (b)

coefficients(reglin) # mostra i coefficienti dell'equazione hc = a + b · rcc

confint(reglin, level=0.95) # calcola gli intervalli di confidenza dell'intercetta e del coefficiente angolare

#

# analisi statistica della adeguatezza della regressione lineare

#

summary(reglin) # mostra un riepilogo dei risultati

t.test(residuals(reglin)) # verifica che la media degli errori non sia significativamente diversa da zero

shapiro.test(residuals(reglin)) # verifica la normalità della distribuzione degli errori

summary(gvlma(reglin)) # test globale per l'assunto di linearità

outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers)

#

# analisi grafica della adeguatezza della regressione lineare

#

windows() # apre una nuova finestra

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

newx = seq(min(var_x), max(var_x), by = 0.01) # valori della x per i quali calcolare l'intervallo di confidenza

conf_interval <- predict(reglin, newdata=data.frame(var_x=newx), interval="confidence", level = 0.95) # calcola gli intervalli di confidenza

plot(var_x, var_y, xlab="Altezza dei padri (pollici)", ylab="Altezza dei figli (pollici)", main="Regressione lineare y = a + b·x") # grafico dei dati

abline(reglin, col="lightblue") # retta di regressione

lines(newx, conf_interval[,2], col="blue", lty=2) # limite di confidenza inferiore

lines(newx, conf_interval[,3], col="blue", lty=2) # limite di confidenza superiore

#

plot(var_y, residuals(reglin), xlab="Altezza dei figli (pollici)", ylab="Altezza osservata - calcolata (pollici)", main="Analisi delle differenza residue") # grafico delle differenza tra altezza osservata e altezza calcolata con l'equazione della retta

#

influencePlot(reglin, fill=FALSE, xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

#

qqPlot(reglin, xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

#

Anche in questo caso, per quanto concerne l'analisi statistica, le principali conclusioni sono quelle tratte con il test globale per l'assunto di linearità effettuato mediante la funzione gvlma(), che però questa volta ci dice che gli assunti che stanno alla base del modello di regressione lineare, a parte la curtosi, non sono per niente soddisfatti:

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS

USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:

Level of Significance =  0.05 

Call:

 gvlma(x = reglin) 

                    Value   p-value                   Decision

Global Stat        25.489 4.011e-05 Assumptions NOT satisfied!

Skewness            8.979 2.731e-03 Assumptions NOT satisfied!

Kurtosis            1.965 1.610e-01    Assumptions acceptable.

Link Function       4.632 3.138e-02 Assumptions NOT satisfied!

Heteroscedasticity  9.913 1.641e-03 Assumptions NOT satisfied!

Il test di Bonferroni non evidenzia dati anomali, ma segnala il dato numero 1 come quello che maggiormente si discosta dai rimanenti:

> outlierTest(reglin) # valore p di Bonferonni per la presenza di dati anomali (outliers) 

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

   rstudent unadjusted p-value Bonferonni p

1 -3.512671         0.00046509      0.43161

Le conclusioni dell'analisi grafica confermano, con i valori delle differenze residue che mostrano una proporzionalità diretta con l'altezza, che i dati in questione violano uno degli assunti del modello di regressione lineare 

e confermano anche la possibile presenza di dati anomali con il grafico delle distanze di Cook e con il grafico dei quantili, dati che sono evidenziati sia direttamente nei grafici, sia nei rispettivi riepiloghi riportati nella Console di R:

> influencePlot(lm(galton$child ~ galton$parent), xlab="t-quantili (il diametro dei cerchi", sub="è proporzionale alla distanza di Cook)", ylab="Residui studentizzati", main="Influenza dei dati") # grafico dell'influenza dei dati sulle conclusioni

       StudRes         Hat       CookD

1   -3.5126706 0.002699826 0.016499436

2   -2.9226403 0.001090010 0.004622768

857  0.4839886 0.008511029 0.001006222

897  2.6611087 0.003740530 0.013207269

898  0.9327550 0.008511029 0.003734739

> #

> qqPlot(lm(galton$child ~ galton$parent), xlab="t-quantili", ylab="Residui studentizzati", main="Quantili vs. residui") # mostra il grafico dei quantili per i residui studentizzati

[1] 1 2

I dati numero 1, 2, 857, 897, 898 forniti dalla funzione influencePlot() stanno di nuovo a indicare i casi che influenzano in modo importante la regressione, casi dei quali sarebbe importante controllare la validità, o che potrebbero essere situati in intervalli di valori per i quali potrebbe essere opportuno acquisire più dati. Anche la funzione qqPlot() fornisce, oltre al grafico, l'indicazione di due punti da controllare, che sono punti 1, 2 già indicati dalla funzione precedente.

Il set di dati galton, oltre a non soddisfare i requisiti fondamentali che si richiedono ai dati per l'applicazione della regressione lineare ordinaria, è anche un esempio di dati nei quali non è chiaro, a dispetto delle conclusioni tratte da Francis Galton di una "regressione verso la media" delle altezze dei figli rispetto a quelle dei padri, quale delle due variabili debba essere considerata la variabile indipendente. Le implicazioni di questo fatto, le conseguenze che esso determina nelle conclusioni tratte dalla regressione lineare, e il calcolo della regressione lineare con modelli alternativi, sono discussi a parte nel post La regressione lineare: assunti e modelli.

Come potete notare entrambi gli script sono stato realizzati in modo da rendere immediato il loro riutilizzo: è sufficiente assegnare alla variabile x (var_x <-) e alla variabile y (var_y <-) i vostri nuovi dati (e personalizzare opportunamente titoli e legende).

Per una guida rapida all'importazione dei dati potete consultare i link:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestione dei dati mancanti

----------

[1] Bossi A, Cortinovis I, Duca PG, Marubini E. Introduzione alla statistica medica. La Nuova Italia Scientifica, Roma, 1994, ISBN 88-430-0284-8. Il modello statistico nella regressione lineare, pp-305-308.

[2] Snedecor GW, Cochran WG. Statistical Methods. The Iowa State University Press, 1980, ISBN 0-8138-1560-6. The matematical model in linear regression, pp.153-157.

[3] Regressione lineare.

https://it.wikipedia.org/wiki/Regressione_lineare

[4] Vedere il post La regressione lineare: assunti e modelli. 

[5]  Vedere il post Regressione lineare semplice non parametrica.

[6] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[7] Digitate help(lm) nella Console di R per la documentazione della funzione lm().

[8] Vedere il post Il set di dati galton.

[9] Informazioni esaurienti sono contenute nei manuali di riferimento dei pacchetti aggiuntivi qui impiegati, che trovate alla pagina Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html