Come impiegare l'operatore %>% e l'argomento FUN

Anche se l'obiettivo di questi post rimane sempre quello di fornire gli strumenti più semplici, quelli utili nella fase iniziale di apprendimento – notoriamente tutta in salita e nella quale quello che conta è "rompere il ghiaccio" con R – occasionalmente può essere interessante affrontare argomenti un poco più tecnici, come i due che vediamo ora.

Il primo riguarda l'operatore pipe indicato con il simbolo %>%: il termine sta per il sostantivo "tubo" ma anche per l'espressione verbale "trasportare con una tubatura", ed è proprio in questo senso che viene impiegato. L'operatore pipe è stato introdotto la prima volta con il pacchetto magrittr [1] che ne riporta la seguente descrizione: 

"[il pacchetto] fornisce un meccanismo per concatenare i comandi con il nuovo operatore di inoltro pipe, %>%. Questo operatore inoltrerà un valore, o il risultato di un'espressione, alla successiva chiamata/espressione di funzione. [l'operatore] fornisce un supporto flessibile per il tipo di espressioni sul lato destro. Per ulteriori informazioni, consultare la descrizione del pacchetto. Per citare René Magritte, 'Ceci n'est pas une pipe'. " [2]

In altre parole: l'operatore %>% fornisce la "tubatura" che consente – in modo completamente trasparente per il programmatore – di inoltrare i dati in uscita da una funzione all'ingresso della funzione che segue. Vediamo questo breve esempio. Abbiamo il set di dati ais [3] 

> ais

     rcc  wcc   hc   hg ferr   bmi   ssf pcBfat    lbm    ht    wt sex   sport

1   3.96  7.5 37.5 12.3   60 20.56 109.1  19.75  63.32 195.9  78.9   f  B_Ball

…...

202 5.38  6.3 46.0 15.7   32 21.07  34.9   6.26  72.00 190.8  76.7   m  Tennis

 

e vogliamo calcolare la media, la deviazione standard e l'errore standard della concentrazione di emoglobina (hg) nel sangue dei 202 atleti australiani ivi inclusi suddividendoli per sport praticato (sport) e per sesso (sex).

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

# COME IMPIEGARE l'operatore %>% (pipe)

#

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

library(dplyr) # carica il pacchetto con l'operatore pipe e altre funzioni

library(plotrix) # carica il pacchetto per il calcolo dell'errore standard

# 

# calcola media, deviazione standard ed errore standard separatamente per sport e per sesso

# 

pipe_hg <- ais %>% group_by(sport, sex) %>% summarise(media=mean(hg), ds=sd(hg), es=std.error(hg)) %>% ungroup()

#

pipe_hg # mostra i risultati

#

Dopo avere installato e caricato i pacchetti necessari, per realizzare la tabella che contiene i risultati desiderati impieghiamo più volte l'operatore pipe (%>%), che trasferisce i dati da una funzione alla successiva in questo modo:

→ i dati (ais) sono trasferiti (%>%) alla funzione group_by();

→ la funzione group_by() – che raggruppa i dati in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex) – trasferisce il risultato di questa operazione (%>%) alla funzione summarise();

→ la funzione summarise() calcola la media mean(), la deviazione standard sd() e l'errore standard std.error() della variabile hg – la concentrazione nel sangue dell'emoglobina (espressa in g/dL, grammi per decilitro di sangue) dei 202 atleti australiani suddivisi per sport e per sesso – e trasferisce il risultato di questa operazione (%>%) alla funzione ungroup();

→ la funzione ungroup() rimuove una serie di valori impiegati provvisoriamente per generare il risultato finale, che viene infine salvato (<-) nella tabella pipe_hg.

Ed ecco la tabella risultante, ottenuta quindi con una sola riga di codice:

> pipe_hg # mostra i risultati

# A tibble: 17 × 5

   sport   sex   media    ds    es

   <fct>   <fct> <dbl> <dbl> <dbl>

 1 B_Ball  f      13.1 0.878 0.243

 2 B_Ball  m      15.1 0.922 0.266

 3 Field   f      14.6 0.682 0.258

 4 Field   m      16.0 0.805 0.232

 5 Gym     f      13.6 0.860 0.430

 6 Netball f      12.8 0.567 0.118

 7 Row     f      14.0 0.740 0.158

 8 Row     m      15.4 0.711 0.184

 9 Swim    f      13.6 0.583 0.194

10 Swim    m      15.5 0.655 0.182

11 T_400m  f      13.8 1.04  0.312

12 T_400m  m      15.3 0.824 0.194

13 T_Sprnt f      14.2 0.556 0.278

14 T_Sprnt m      16.2 1.49  0.450

15 Tennis  f      13.5 1.10  0.414

16 Tennis  m      15.6 1.48  0.741

17 W_Polo  m      15.5 0.718 0.174

Il secondo riguarda l'argomento FUN che consente di specificare una o più funzioni e di applicarle in sequenza a una variabile. Impieghiamo gli stessi dati dell'esempio precedente allo stesso scopo: calcolare la media, la deviazione standard e l'errore standard della concentrazione di emoglobina nel sangue di 202 atleti australiani suddividendoli per sport praticato e per sesso. 

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

# COME IMPIEGARE l'argomento FUN

#

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

library(plotrix) # carica il pacchetto per il calcolo dell'errore standard

#

# calcola media, deviazione standard ed errore standard separatamente per sport e per sesso

#

FUN_hg <- aggregate(hg~sport+sex, data=ais, FUN=function(x) c(media=mean(x), ds=sd(x), es=std.error(x))) 

#

FUN_hg # mostra i risultati 

#

Dopo avere installato e caricato i pacchetti necessari, la tabella che contiene i risultati desiderati viene realizzata con la funzione aggregate() che nell'ordine:

→ nel nostro set di dati (data=ais) raggruppa (~) i valori di concentrazione dell'emoglobina (hg) in sottoinsiemi per sport praticato (sport) e per sesso dell'atleta (sex);

→ con l'argomento FUN e con il valore =function(x) specifica che sui dati devono essere effettuati i calcoli di seguito elencati;

→ con la funzione c() elenca le funzioni da impiegare per i calcoli;

→ con le funzioni mean(), sd() e std.error() effettua i calcoli di media, deviazione standard ed errore standard sui dati raggruppati per sport e per sesso, salvando i risultati nelle variabili media, ds ed es. 

I risultati sono infine riportati (<-) nella tabella FUN_hg e, a parte l'ordinamento e gli arrotondamenti, sono identici a quelli ottenuti con lo script precedente, impiegando anche in questo caso una sola riga di codice.

> FUN_hg # mostra i risultati 

     sport sex   hg.media      hg.ds      hg.es

1   B_Ball   f 13.1307692  0.8778616  0.2434750

2    Field   f 14.6285714  0.6824326  0.2579353

3      Gym   f 13.6000000  0.8602325  0.4301163

4  Netball   f 12.8173913  0.5670114  0.1182301

5      Row   f 14.0318182  0.7396179  0.1576871

6     Swim   f 13.5666667  0.5830952  0.1943651

7   T_400m   f 13.7727273  1.0354621  0.3122036

8  T_Sprnt   f 14.1750000  0.5560276  0.2780138

9   Tennis   f 13.5285714  1.0965313  0.4144499

10  B_Ball   m 15.1416667  0.9219134  0.2661335

11   Field   m 16.0250000  0.8046738  0.2322893

12     Row   m 15.3866667  0.7110020  0.1835799

13    Swim   m 15.5076923  0.6550592  0.1816807

14  T_400m   m 15.3055556  0.8235306  0.1941080

15 T_Sprnt   m 16.1909091  1.4936228  0.4503442

16  Tennis   m 15.6500000  1.4821156  0.7410578

17  W_Polo   m 15.5176471  0.7178399  0.1741017

Conclusione: nel corso dell'elaborazione dei dati può essere utile tenere a portata di mano qualche esempio di codice, come questi due, che ricorda come impiegare l'operatore %>% e l'argomento FUN. Non trovate molto altro sul web, ma potete provare ricercando ad esempio "R %>% operator examples" o "R FUN function examples".

-----------

[1] Vedere il Reference manual del pacchetto: Package ‘magrittr’.

https://cran.r-project.org/web/packages/magrittr/magrittr.pdf

[2] "Provides a mechanism for chaining commands with a new forward-pipe operator, %>%. This operator will forward a value, or the result of an expression, into the next function call/expression. There is flexible support for the type of right-hand side expressions. For more information, see package vignette. To quote Rene Magritte, 'Ceci n'est pas un pipe.' ".

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

Valutare l'omogeneità tra varianze

L'analisi della varianza (ANOVA) è basata "... su due importanti ipotesi: (a) la normalità delle distribuzioni delle osservazioni … e (b) la costanza delle varianze nei diversi gruppi" [1]. 

Per valutare la normalità (gaussianità) delle distribuzioni vi è solo l'imbarazzo della scelta tra i molti test disponibili, molti dei quali li trovate illustrati anche in questo sito [2]. 

Poiché invece qualche difficoltà può sorgere in merito al metodo per valutare la omogeneità delle varianze, qui illustro come per questo si può impiegare il test di Levene.

I dati sono quelli già impiegati nella ANOVA a un fattore [3], copiate le sedici righe riportate qui sotto 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).

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

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

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

# Test di Levene per l'omogeneità tra varianze

#

library(car) # carica il pacchetto necessario per eseguire il test

#

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

#

leveneTest(produzione~macchina, data=mydata) # verifica se le varianze risultano omogenee

#

Dopo avere caricato il pacchetto car che include la funzione leveneTest() sono importati i dati nell'oggetto mydata. 

Il senso dello script risiede tutto nella terza riga di codice: prima di eseguire la ANOVA vogliamo verificare in via preliminare se le varianze dei cinque dati di produzione rilevati per ciascuna delle tre macchine i1, i2 e i3 sono tra loro omogenee applicando il test di Levene sui dati (mydata) con i valori di produzione aggregati per macchina (produzione~macchina). Questo è quanto ci compare:

> leveneTest(produzione~macchina, data=mydata) # verifica se le varianze risultano omogenee

Levene's Test for Homogeneity of Variance (center = median)

      Df F value Pr(>F)

group  2  0.4238  0.664

      12               

Messaggio di avvertimento:

In leveneTest.default(y = y, group = group, ...) : group coerced to factor.

Perché le tre varianze calcolate per le produzioni delle tre macchine possano essere considerate diverse dovremmo avere p<0.05, ma qui abbiamo un p=0.664 e pertanto concludiamo che le varianze delle produzioni ottenute con le tre diverse macchine non differiscono tra loro in modo significativo. Quindi l'omogeneità delle varianze – uno dei due prerequisiti della ANOVA – risulta rispettato.

Da notare che il Messaggio di avvertimento appare in quanto la variabile di raggruppamento utilizzata per calcolare il test di Levene deve essere una variabile "fattore". Se la variabile utilizzata non lo è, ma è comunque la variabile qualitativa con la quale vogliamo raggruppare i dati, come "macchina" nel nostro esempio, la funzione leveneTest() ci viene in aiuto in quanto trasforma automaticamente la nostra variabile in una variabile "fattore", pertanto il messaggio può essere ignorato.

Conclusione: i risultati dell'ANOVA devono essere validati accertando che siano rispettati due requisiti:

(a) i dati devono essere distribuiti in modo gaussiano, cosa che può essere verificata impiegando qualcuno dei numerosi test di normalità (gaussianità) disponibili;

(b) i dati devono avere varianze omogenee (varianze non significativamente diverse), cosa che può essere verificata impiegando il test di Levene qui illustrato.

----------

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

[2] I metodi per l'analisi della normalità (gaussianità) di una distribuzione sono riportati in:

- Test di normalità (gaussianità);

- Valutare asimmetria e curtosi;

- Tabulare una serie di test di normalità (gaussianità).

[3] Vedere il post Analisi della varianza a un fattore.

Analisi delle distribuzioni dei dati nei gruppi/cluster

Se trovate utile l'impiego delle tecniche di clustering (o analisi dei gruppi) [1] come metodo per l'analisi esplorativa dei dati – ma anche se le affrontate per la prima volta – potrebbe interessarvi questo esempio.

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

# ANALIZZARE LE DISTRIBUZIONI DEI DATI NEI CLUSTER con ggpairs  

#  

library(cluster) # carica il pacchetto per il clustering

library(factoextra) # carica un pacchetto addizionale per il clustering

library(GGally) # carica il pacchetto che estende le funzioni di ggplot2

#

# preparazione dei dati

#

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

newdata <- subset(mydata, select=c("gr", "hb", "mcv", "hb_a2", "ferro")) # seleziona le variabili per il clustering

#

# clustering non gerarchico con il metodo di Rousseew (k-medoids)

#

myclust <- pam(newdata, 3, metric=c("euclidean"), stand=TRUE) # genera un oggetto 'pam' che contiene i risultati del clustering

fviz_cluster(myclust, myclust$cluster, labelsize=7, main = "Grafico dei cluster - metodo di Rousseew (k-medoids)") # grafico dei cluster per le prime due componenti principali

#

Vediamo di cosa si tratta. Abbiamo nel file ema.csv una raccolta di 643 casi dei quali abbiamo raccolto i risultati di tre analisi del sangue richieste a scopo diagnostico:

→ esame emocromocitometrico – del quale ci interessano in particolare i valori di: (i) gr la concentrazione degli eritrociti nel sangue (globuli rossi) espressa in migliaia di miliardi per litro di sangue 10¹²/L o in milioni per microlitro di sangue 10⁶/µL (le due espressioni sono numericamente equivalenti); (ii) hb la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue g/dL; (iii) mcv il volume medio dei globuli rossi espresso in milionesimi di miliardesimo di litro (10-¹⁵) o femtolitri (fL);

→ emoglobina A2 – hb_A2 la variante fisiologica dell'emoglobina presente nell'adulto espressa come percentuale dell'emoglobina totale;

→ ferro – ferro la concentrazione del ferro espressa in microgrammi per decilitro di siero µg/dL. 

Sappiamo che i 643 casi comprendono:

→ soggetti sani, impiegati come controlli, con i risultati delle tre analisi all'interno dei valori fisiologici (valori "normali");

→ soggetti con alterazioni dell'esame emocromocitometrico e con carenza di ferro (confermata dalla diminuita concentrazione del ferro nel siero);

→  soggetti con alterazioni dell'esame emocromocitometrico e portatori del gene della beta-talassemia (confermato dall'aumento della percentuale di emoglobina A2);

→  soggetti con alterazioni dell'esame emocromocitometrico ma non classificabili ulteriormente con i dati disponibili.

Siamo interessati a una analisi esplorativa dei dati e, seguendo lo script, procediamo come segue.

1) Preparazione dei dati. 

Installiamo i tre pacchetti che ci servono, cluster, factoextra e GGally, quindi ci procuriamo i dati impiegati nell'esempio (per il file di dati trovate link e istruzioni alla pagina di Download):

→ effettuando il download del file di dati ema.csv

→ salvando il file nella cartella C:\Rdati\

Con la funzione read.table() importiamo i dati nell'oggetto mydata. Quindi con la funzione subset() riportiamo nell'oggetto newdata le sole variabili che ci interessano, cioè i risultati delle analisi indicate sopra (643 record), in cinque colonne:

> newdata

      gr   hb  mcv hb_a2 ferro

1   4.90 13.3 82.8   1.8   106

2   4.66 10.8 73.6   2.6   148

3   5.43 11.5 66.5   4.8   104

4   5.41 10.8 73.4   2.5    74

..........

642 5.28 11.0 65.4   1.9    11

643 4.03  7.2 55.0   1.9    21

2) Clustering non gerarchico (analisi dei gruppi). 

A questo punto eseguiamo il clustering non gerarchico, lo facciamo impiegando il metodo di Rousseew. La prima funzione utilizzata è pam() che deriva il suo nome da "Partitioning Around Medoids", impiega l'algoritmo di clustering dei dati  k-medoids (un versione robusta dell'algoritmo k-means) e prevede come argomenti:

→ i dati newdata da impiegare;

→ il numero di cluster, in questo caso 3, in cui classificare i dati;

→ la metrica da impiegare che è "euclidean" (in alternativa può essere "manhattan");

→ stand=TRUE che impone la standardizzazione dei dati, che viene eseguita automaticamente dalla funzione.

Con la successiva funzione fviz_cluster() è generato il grafico dei cluster (myclust$cluster) ottenuti con l'algoritmo k-medoids e contenuti nell'oggetto myclust.

Date le premesse riportate all'inizio abbiamo impostato 3 come numero dei gruppi (cluster) in cui classificare i dati allo scopo di separare (i) sani, (ii) carenza di ferro, (iii) beta-talassemia eterozigote. Consapevoli del fatto che data la casistica, la quantità e il tipo delle informazioni di cui disponiamo, le altre possibili e meno frequenti cause di alterazione dell'esame emocromocitometrico (e ve ne sono molte) rimarranno inestricabilmente sovrapposte tra loro e/o con qualcuno di questi tre gruppi (a causa dei vincoli imposti dalla informazione disponibile che è, come sempre, limitata).

3) Analisi delle distribuzioni dei dati nei cluster.

Eccoci arrivati al punto. Non ci basta avere, per ciascun cluster, l'indicazione dei soggetti (numerati in ordine progressivo in mydata) appartenenti a ciascun cluster. Vogliamo avere anche l'indicazione della corrispondenza tra i cluster e le distribuzioni dei risultati delle grandezze misurate. Lo facciamo in una nuova finestra grafica aperta con windows() associando a ciascun caso (record) della tabella newdata con la funzione cbind() il cluster di appartenenza myclust$clustering estratto dall'oggetto 'pam' myclust che contiene i risultati del clustering e che viene riportato aggiungendo una sesta colonna alla tabella.

#

# analisi delle distribuzioni dei dati nei cluster

#

windows() # apre e inizializza una nuova finestra grafica

newclust <- cbind(newdata, myclust$clustering) # riporta per ciascun caso il cluster di appartenenza

ggpairs(newclust, columns=1:5, aes(color=as.factor(myclust$clustering), alpha=0.5), title="Scatterplot e kernel density plot dei dati dei cluster") # scatterplot e kernel density plot dei dati nei cluster

# 

Ora abbiamo la tabella che possiamo impiegare per generare il grafico della distribuzione dei dati dei cluster con la funzione ggpairs() che impiega gli argomenti:

→ newclust per indicare la tabella che contiene i dati;

→ 1:5 per indicare che i dati da rappresentare sono contenuti nelle colonne da 1 a 5 (la sesta contiene il cluster di appartenenza);

→ as.factor() in quanto le variabili numeriche 1, 2 e 3 della colonna 6 che indicano il cluster di appartenenza devono essere trasformate in variabili fattore;

→ alpha=0.5 per rappresentare i colori con una trasparenza del 50%.

Da notare che al bisogno potete anche esportare la tabella newclust, che nelle prime cinque colonne contiene i dati e nella sesta il cluster di appartenenza come qui indicato 

> newclust

      gr   hb  mcv hb_a2 ferro myclust$clustering

1   4.90 13.3 82.8   1.8   106                  1

2   4.66 10.8 73.6   2.6   148                  1

3   5.43 11.5 66.5   4.8   104                  2

4   5.41 10.8 73.4   2.5    74                  3

..........

642 5.28 11.0 65.4   1.9    11                  3

643 4.03  7.2 55.0   1.9    21                  3

e salvarla sotto forma di file [2].

Ed ecco il grafico risultante.

La conclusione? 

1) l'analisi dei gruppi fornisce una sintesi nella quale le grandezze soggiacenti non sono più riconoscibili singolarmente in quanto le componenti principali sono una combinazione lineare di tutte le variabili originali (nel primo dei grafici la prima componente principale in ascisse è la variabile formata dalla combinazione lineare di tutte le variabili originali che spiega la maggior quantità di varianza, la seconda componente principale in ordinate è la variabile che spiega la maggior quantità di varianza in ciò che rimane una volta rimosso l'effetto della prima componente principale);

2) la distribuzione dei dati all'interno dei cluster fornisce una rappresentazione di dettaglio che consente di ricollegare questi ultimi con la distribuzione originaria delle variabili la cui combinazione ha dato loro origine, e consente di riconoscere come:

→ il cluster 1 (rosso) contiene i soggetti con emoglobina A2 e ferro normali;

→ il cluster 2 (verde) contiene i soggetti con emoglobina A2 aumentata; 

→ il cluster 3 (azzurro) contiene i soggetti con ferro diminuito.

Pur con distribuzioni che in parte si sovrappongono – un fatto determinato dalla incompletezza della informazione disponibile rispetto alla complessità dei casi – l'esempio illustra come:

→ il clustering fornisce una vista sintetica dell'informazione globale contenuta nei dati;

→ possiamo aggiungere all'analisi dei cluster una vista di dettaglio delle distribuzioni dei singoli dati all'interno dei cluster;

→ le due viste degli stessi dati si integrano l'una con l'altra. 

----------

[1] Fate click su clustering nelle Parole chiave.

[2] Vedere il post Esportazione dei dati in un file .csv.

Analisi della varianza a due fattori con replicati

L'analisi della varianza o ANOVA – denominazione tradizionalmente impiegata per indicare una serie di tecniche che, a dispetto del nome, consentono di effettuare un confronto globale tra più medie – è stata presentata e illustrata in tutti i suoi aspetti in due post ai quali si rimanda [1].

Qui impieghiamo lo stesso disegno sperimentale riportato per l'analisi della varianza a due fattori e cioè

ma questa volta ciascuna delle misure di produzione per macchina/operatore è stata ripetuta quattro volte al fine di migliorarne la stima. Così ad esempio per macchina i=1 e operatore j=2 nel file di dati che impieghiamo oltre al valore singolo 45.7 riportato qui sopra nella tabella originaria abbiamo in aggiunta le misure 46.9, 49.1 e 45.2 (questi dati sono stati evidenziati in colore nei dati riportati qui sotto) e quindi in totale quattro misure replicate. Questo è il disegno sperimentale tipico della analisi della varianza a due fattori con replicati.

Per proseguire è necessario:

→ effettuare il download del file di dati anova2r.csv

→ salvare il file nella cartella C:\Rdati\

Per il file di dati trovate link e modalità di download alla pagina Dati, ma potete anche semplicemente copiare i dati riportati qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarli in C:\Rdati\ in un file di testo denominato anova2r.csv (assicuratevi che il file sia effettivamente salvato con l'estensione .csv).

macchina;operatore;replicato;produzione

i1;j1;1;57.2

i1;j1;2;56.7

i1;j1;3;55.9

i1;j1;4;58.3

i1;j2;1;45.7

i1;j2;2;46.9

i1;j2;3;49.1

i1;j2;4;45.4

i1;j3;1;47.3

i1;j3;2;49.4

i1;j3;3;48.1

i1;j3;4;47.5

i1;j4;1;54.9

i1;j4;2;52.6

i1;j4;3;53.5

i1;j4;4;54.1

i1;j5;1;35.9

i1;j5;2;37.2

i1;j5;3;38.7

i1;j5;4;36.4

i2;j1;1;64.5

i2;j1;2;61.1

i2;j1;3;63.9

i2;j1;4;64.7

i2;j2;1;52.4

i2;j2;2;55.3

i2;j2;3;57.2

i2;j2;4;53.9

i2;j3;1;54.1

i2;j3;2;56.2

i2;j3;3;52.3

i2;j3;4;53.8

i2;j4;1;57.5

i2;j4;2;55.9

i2;j4;3;56.3

i2;j4;4;58.8

i2;j5;1;52.5

i2;j5;2;50.3

i2;j5;3;51.4

i2;j5;4;53.3

i3;j1;1;56.6

i3;j1;2;59.7

i3;j1;3;56.3

i3;j1;4;56.5

i3;j2;1;51.6

i3;j2;2;50.3

i3;j2;3;47.8

i3;j2;4;52.6

i3;j3;1;49.5

i3;j3;2;52.7

i3;j3;3;49.1

i3;j3;4;47.7

i3;j4;1;56.9

i3;j4;2;58.5

i3;j4;3;57.1

i3;j4;4;55.2

i3;j5;1;44.2

i3;j5;2;46.7

i3;j5;3;43.9

i3;j5;4;42.6

Per l'analisi della varianza anche in questo caso viene impiegata la funzione aov() impiegata per le altre [1] specificando questa volta che la variabilità deve essere decomposta in variabilità dovuta alla macchina, dovuta all'operatore e dovuta ai replicati eseguiti. 

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

# ANOVA analisi della varianza a due fattori con replicati

#

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

#

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

#

summary(anova2r) # mostra i risultati

#

Il risultato della ANOVA

> summary(anova2r) # mostra i risultati

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

macchina     2  602.8   301.4  54.792 1.58e-13 ***

operatore    4 1552.2   388.1  70.544  < 2e-16 ***

replicato    1    0.3     0.3   0.048    0.827    

Residuals   52  286.0     5.5                     

---

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

ci dice che le differenze di produzione tra le tre macchine e le differenze di produzione tra i cinque operatori sono significative, mentre le misure di produzione replicate all'interno della stessa coppia macchina/operatore non differiscono significativamente (p=0.827).

Possiamo avere anche un riscontro grafico dei dati, copiate e incollate nella Console di R queste righe e premete ↵ Invio.

# traccia boxplot per macchina e operatore

#

boxplot(produzione~macchina+operatore, data=mydata) 

#

(per mostrare sulle ascisse i nomi di tutte le variabili, dopo avere generato il grafico ho "afferrato" con il mouse il lato sinistro della finestra per allargarla)

A questo punto manca però la verifica che i dati siano conformi ai due assunti alla base dell'ANOVA, che abbiano cioè:

→ una distribuzione normale (gaussiana);

→ varianze omogenee;

due esigenze che abbiamo già illustrato nei precedenti post sulla ANOVA [1].

Per valutare la normalità (gaussianità) dei dati impieghiamo il test di normalità di Shapiro-Wilk, copiate queste ultime righe, incollatele nella Console di R e premete ↵ Invio:

# ANOVA verifica della normalità delle distribuzioni

#

mac <- split(mydata$produzione, mydata$macchina) # estrae e separa i dati di produzione per macchina

#

sapply(mac, shapiro.test) # verifica se i dati delle macchine sono distribuiti in modo gaussiano

#

La funzione split() estrae dall'oggetto mydata i dati di produzione (mydata$produzione) per ciascuna macchina (mydata$macchina) e li riporta, separatamente, nell'oggetto mac. Se digitate

mac 

potete vedere come sono organizzati i dati estratti per le tre macchine.

> mac

$i1

 [1] 57.2 56.7 55.9 58.3 45.7 46.9 49.1 45.4 47.3 49.4 48.1 47.5 54.9 52.6 53.5

[16] 54.1 35.9 37.2 38.7 36.4

$i2

 [1] 64.5 61.1 63.9 64.7 52.4 55.3 57.2 53.9 54.1 56.2 52.3 53.8 57.5 55.9 56.3

[16] 58.8 52.5 50.3 51.4 53.3

$i3

 [1] 56.6 59.7 56.3 56.5 51.6 50.3 47.8 52.6 49.5 52.7 49.1 47.7 56.9 58.5 57.1

[16] 55.2 44.2 46.7 43.9 42.6

A questo punto con la funzione sapply() il test di Shapiro-Wilk viene applicato automaticamente ai dati delle tre macchine.

> sapply(mac, shapiro.test) # verifica se i dati delle macchine sono distribuiti in modo gaussiano

          i1                            i2                           

statistic 0.9202174                     0.9052156                    

p.value   0.1000444                     0.05170547                   

method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"

data.name "X[[i]]"                      "X[[i]]"                     

          i3                           

statistic 0.9449654                    

p.value   0.2970447                    

method    "Shapiro-Wilk normality test"

data.name "X[[i]]"   

I valori di p ci dicono che le distribuzioni possono essere considerate gaussiane, quindi questo assunto dell'ANOVA è rispettato. 

Ora valutiamo l'omogeneità delle varianze impiegando il test di Levene, copiate queste righe, incollatele nella Console di R e premete ↵ Invio:

# Test di Levene per l'omogeneità tra varianze

#

library(car) # carica il pacchetto necessario per eseguire il test

#

leveneTest(produzione~macchina, data=mydata) # verifica se le varianze risultano omogenee

#

Dopo avere caricato il pacchetto car il test di Levene viene eseguito impiegando la funzione leveneTest() sui dati di produzione per ciascuna delle tre macchine (produzione~macchina). Il messaggio di avvertimento può essere ignorato in quanto ci dice semplicemente che la funzione ha trasformato automaticamente la nostra variabile in una variabile "fattore".

> leveneTest(produzione~macchina, data=mydata) # verifica se le varianze risultano omogenee

Levene's Test for Homogeneity of Variance (center = median)

      Df F value  Pr(>F)  

group  2  2.4955 0.09142 .

      57                  

---

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

Messaggio di avvertimento:

In leveneTest.default(y = y, group = group, ...) : group coerced to factor.

Il valor di p (0.09142) ci dice che le varianze sono omogenee, quindi anche questo assunto è rispettato.

La conclusione: stabilito che per i dati di produzione per macchina i due requisiti del modello statistico alla base dell'ANOVA sono rispettati, possiamo accettare il risultato della analisi della varianza a due fattori con replicati che ci dice che esiste una differenza statisticamente significativa tra le medie di produzione delle tre macchine (p=1.58e-13).

Sono lasciati come esercizio l'applicazione del test di normalità e del test di omogeneità tra varianze ai dati di produzione per operatore e le relative conclusioni, mentre per questo specifico disegno della ANOVA valgono, identiche, le considerazioni generali riportate nei due post precedenti, ai quali si rimanda [1]

----------

[1] Vedere il post Analisi della varianza a un fattore e il post Analisi della varianza a due fattori.