Visualizzazione post con etichetta split(). Mostra tutti i post
Visualizzazione post con etichetta split(). Mostra tutti i post

sabato 21 dicembre 2024

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 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]


----------


venerdì 22 novembre 2024

Estrazione dei dati per fattore

Il set di dati ais contenuto nel pacchetto DAAG – accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [1] – rappresenta un esempio paradigmatico di dati estratti da un database e sui quali si potrebbe voler eseguire una qualche analisi statistica o grafica, lo vediamo e poi lo impieghiamo per illustrare un esempio di estrazione dei dati per fattore.

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

# ESTRAZIONE DEI DATI PER FATTORE
#
library(DAAG) # carica il pacchetto DAAG che include il set di dati ais
head(ais, 3) # mostra le prime tre righe del set di dati ais
tail(ais,3) # mostra le ultime tre righe del set di dati ais
str(ais) # mostra la struttura dei dati 
#

Dopo avere caricato il pacchetto, con head(ais, 3) vediamo le prime tre righe di ais

> head(ais, 3) # mostra le prime tre righe del set di dati 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
2 4.41 8.3 38.2 12.7   68 20.67 102.8  21.30 58.55 189.7 74.4   f B_Ball
3 4.14 5.0 36.4 11.6   21 21.86 104.6  19.88 55.36 177.8 69.1   f B_Ball

con tail(ais, 3) vediamo le ultime tre righe 

> tail(ais,3) # mostra le ultime tre righe del set di dati ais
     rcc wcc   hc   hg ferr   bmi  ssf pcBfat lbm    ht   wt sex  sport
200 5.03 6.4 42.7 14.3  122 22.01 47.6   8.51  68 183.1 73.8   m Tennis
201 4.97 8.8 43.0 14.9  233 22.34 60.4  11.50  63 178.4 71.1   m Tennis
202 5.38 6.3 46.0 15.7   32 21.07 34.9   6.26  72 190.8 76.7   m Tennis

e con str(ais) vediamo la struttura dei dati contenuti in ais

> str(ais) # mostra la struttura dei dati 
'data.frame':   202 obs. of  13 variables:
 $ rcc   : num  3.96 4.41 4.14 4.11 4.45 4.1 4.31 4.42 4.3 4.51 ...
 $ wcc   : num  7.5 8.3 5 5.3 6.8 4.4 5.3 5.7 8.9 4.4 ...
 $ hc    : num  37.5 38.2 36.4 37.3 41.5 37.4 39.6 39.9 41.1 41.6 ...
 $ hg    : num  12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...
 $ ferr  : num  60 68 21 69 29 42 73 44 41 44 ...
 $ bmi   : num  20.6 20.7 21.9 21.9 19 ...
 $ ssf   : num  109.1 102.8 104.6 126.4 80.3 ...
 $ pcBfat: num  19.8 21.3 19.9 23.7 17.6 ...
 $ lbm   : num  63.3 58.5 55.4 57.2 53.2 ...
 $ ht    : num  196 190 178 185 185 ...
 $ wt    : num  78.9 74.4 69.1 74.9 64.6 63.7 75.2 62.3 66.5 62.9 ...
 $ sex   : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ...
 $ sport : Factor w/ 10 levels "B_Ball","Field",..: 1 1 1 1 1 1 1 1 1 1 ...

Tutto molto chiaro e molto semplice: abbiamo 202 righe di dati con 10 variabili numeriche (num risultato di misure di dati ematologici e di dati biometrici effettuate in altrettanti atleti australiani [1] – e due variabili "fattori" (Factor) cioè variabili qualitative che sono state impiegate per classificare gli atleti in base al sesso (sex) in due classi (levels) e in base allo sport praticato (sport) in dieci classi. 

Eccoci quindi al punto. Ora che abbiamo l'oggetto sul quale effettuare una analisi statistica dei dati vogliamo fare quello che in genere si fa prima di procedere ulteriormente: capire se i dati sono distribuiti in modo normale (gaussiano), onde applicare i test appropriati. Ma vogliamo anche evitare di ripetere la procedura dieci volte, richiamando le variabili una per una con altrettante righe di codice.

La soluzione pratica che trovo più semplice e che vediamo ora prevede l'impiego della funzione split() che applichiamo alla variabile emoglobina (hg) della tabella ais (ais$hg) per estrarre e separare i dati in base al fattore sport (ais$sport).

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

#
# estrae e separa i dati in base al fattore sport
#
dati_sport <- split(ais$hg, ais$sport) 
#

Per capire cosa è accaduto vediamo l'oggetto dati_sport che abbiamo creato. Copiate e incollate nella Console di R questo codice e premete ↵ Invio.

#
# mostra i dati separati per i dieci valori riscontrati nel fattore sport
#
str(dati_sport) 
#

Ed ecco la struttura del nuovo oggetto, che contiene i valori di concentrazione dell'emoglobina, ma ora separati per i dieci sport praticati.

> str(dati_sport) # per i dieci valori del fattore sport
List of 10
 $ B_Ball : num [1:25] 12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...
 $ Field  : num [1:19] 13.3 14.7 15.3 14.3 14.6 15 15.2 15.4 16.7 15.5 ...
 $ Gym    : num [1:4] 13.4 14 12.5 14.5
 $ Netball: num [1:23] 13.6 12.7 12.3 12.3 12.8 11.8 12.7 12.4 12.4 14.1 ...
 $ Row    : num [1:37] 13.9 14.7 13.3 12.9 14.7 14.3 14.5 13 12.5 14.5 ...
 $ Swim   : num [1:22] 13.8 13.3 12.9 12.7 14.4 14 13.9 14 13.1 15.9 ...
 $ T_400m : num [1:29] 15.9 12.4 14.7 13.9 13 13.8 13.2 12.6 13.5 13.7 ...
 $ T_Sprnt: num [1:15] 13.4 14.4 14.7 14.2 15.6 15.2 15.5 15.5 14.9 19.2 ...
 $ Tennis : num [1:11] 12 13.9 13.5 12.1 14.8 14.5 13.9 17.7 14.3 14.9 ...
 $ W_Polo : num [1:17] 14.4 15.2 15 16.2 15.6 16.2 14.8 15.8 14.4 15.8 ...

Il lavoro preliminare è stato oltremodo rapido, ma quel che ci interessa mostrare è che ora diventa rapidissimo applicare l'analisi a tutte le variabili estratte, basta una sola riga di codice, copiate e incollate nella Console di R questo codice e premete ↵ Invio.

#
# effettua il test di normalità (gaussianità) di Shapiro sui dati dei dieci sport
#
sapply(dati_sport, shapiro.test) 
#

Come si vede dopo avere opportunamente organizzato i dati la funzione sapply() applica il test di Shapiro (shapiro.test) per la normalità (gaussianità) automaticamente a tutte le variabili contenute nell'oggetto dati_sport, cioè alle dieci variabili ricavate dalla scomposizione dei dati in base al fattore sport.

> sapply(dati_sport, shapiro.test) # test di Shapiro sui dati dei dieci sport
          B_Ball                        Field                        
statistic 0.9583881                     0.9740128                    
p.value   0.3833057                     0.8528473                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]"                     
          Gym                           Netball                      
statistic 0.9791217                     0.95727                      
p.value   0.896848                      0.4105184                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]"                     
          Row                           Swim                         
statistic 0.9782721                     0.9457531                    
p.value   0.6714365                     0.2595768                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]"                     
          T_400m                        T_Sprnt                      
statistic 0.9780738                     0.8952224                    
p.value   0.7872988                     0.08047871                   
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]"                     
          Tennis                        W_Polo                       
statistic 0.9399817                     0.9481069                    
p.value   0.520283                      0.4272774                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]"                

Gli stessi dati possono ora essere impiegati anche per generare i grafici a scatola con i baffi delle dieci distribuzioni, copiate e incollate nella Console di R questo codice e premete ↵ Invio.

#
# boxplot dei dati dei dieci sport
#
boxplot(dati_sport) 
#

Ed ecco il grafico, in una versione base ma adeguata (per vedere i nomi di tutte le variabili, dopo avere generato il grafico potete "afferrare" con il mouse il lato sinistro della finestra e allargarla).


L'estrazione e la separazione dei dati in base al fattore sex viene lasciata come esercizio.

La conclusione? Quando si inizia a lavorare con R molte volte ci si arena su problemi banali per i quali esistono soluzioni molto semplici delle quali non si è a conoscenza data la vastità del linguaggio. Un esempio ne è l'estrazione dei dati per "fattore" da un dataset mediante la funzione split(), che consente con estrema semplicità, indicando la variabile da estrarre e il fattore in base al quale separare i dati, di preparare i dati organizzandoli in un formato che può essere molto utile nella successiva analisi esplorativa dei dati stessi.


----------

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