Capitolo 9 Modelli ANOVA con fattori di blocco

Nel capitolo 7 abbiamo visto come è possibile costruire modelli nei quali abbiamo una sola variabile dipendente quantitativa e una variabile indipendente categorica. Abbiamo anche visto che questo tipo di modelli lineari sono normalmente noti come modelli ‘ANOVA’ ad una via e presuppongono che le unità sperimentali siano totalmente indipendenti tra loro. Per questo motivo, essi possono essere utilizzati solo con disegni sperimentali a randomizzazione completa, dove non esistano raggruppamenti di alcun tipo, escluso quello dettato dal trattamento in studio.

Ovviamente è possibile definire modelli analoghi, ma più complessi, nei quali introdurre anche l’effetto di eventuali blocking factors. Vediamo ora alcuni esempi

9.1 Caso-studio: confronto tra erbicidi in campo

Abbiamo una prova di confronto tra erbicidi in mais, con 13 formulati, due testimoni inerbiti (che, per comodità, considereremo due trattamenti diversi) e un testimone scerbato. Date le dimensioni della prova, è lecito ipotizzare che, pur scegliendo un appezzamento il più omogeneo possibile, ci potrebbero essere differenze di infestazione tra un punto e l’altro del campo, con un presumibile gradiente procedendo dai lati (vicino alle fosse) verso il centro. In questa situazione, se l’esperimento fosse disegnato a randomizzazione completa, le differenze di infestazione tra una parte e l’altra del campo sarebbero trascurate e finirebbero per incrementare l’errore sperimentale, diminuendo l’efficienza dell’esperimento.

Si impiega quindi un disegno a blocchi randomizzati con quattro repliche. Ricordiamo che, con questo disegno, il campo è suddiviso in tante sezioni (dette blocchi) quante sono le repliche (quattro), perpendicolarmente al gradiente di infestazione trasversale. In questo modo, l’ambiente è relativamente omogeneo all’interno di ciascun blocco, nel quale viene collocata una replica per trattamento.

Per questa lezione è necessario caricare il package ‘aomisc’, che deve essere preventivamente installato, facendo riferimento, se necessario, al capitolo introduttivo. Il dataset dei risultati (‘rimsulfuron’) è, appunto, contenuto in questo package.

library(aomisc)
data(rimsulfuron)

Più sotto, riportiamo l’output di R relativo ai dati tabulati, con le relative medie di riga (trattamento) e di colonna (blocco)

##                                                 1       2       3       4  Medie
## Alachlor + terbuthylazine                  12.060  49.580  41.340  16.370 29.838
## Hand-Weeded                                77.580  92.080  86.590  99.630 88.970
## Metolachlor + terbuthylazine (pre)         51.770  52.100  49.460  34.670 47.000
## Pendimethalin (post) + rimsuulfuron (post) 94.820  87.720 102.050 101.940 96.632
## Pendimethalin (pre) + rimsulfuron (post)   65.510  88.720  95.520  82.390 83.035
## Rimsulfuron (40)                           85.910  91.090 111.420  93.150 95.392
## Rimsulfuron (45)                           93.030 105.000  89.190  79.040 91.565
## Rimsulfuron (50)                           86.930 105.820 110.020  89.100 97.968
## Rimsulfuron (50+30 split)                  71.360  77.570 115.910  92.160 89.250
## Rimsulfuron (60)                           52.990 102.860 100.620  97.040 88.377
## Rimsulfuron + Atred                        94.110  89.860 104.340  99.630 96.985
## Rimsulfuron + hoeing                       73.220  86.060 118.010  98.320 93.903
## Rimsulfuron + thyfensulfuron               75.280  82.590  94.960  85.850 84.670
## Thifensulfuron                             78.470  42.320  62.520  24.340 51.913
## Unweeded 1                                 10.880  31.770  23.920  20.850 21.855
## Unweeded 2                                 27.580  51.550  25.130  38.610 35.718
## Medie                                      65.719  77.293  83.188  72.068 74.567

9.2 Definizione di un modello lineare

La produzione di ogni unità sperimentale (parcella) è condizionata da più di un effetto:

  1. il diserbante con cui essa è stata trattata;
  2. il blocco di cui essa fa parte;
  3. ogni altro effetto non conoscibile e puramente casuale.

Il modello è quindi:

\[ Y_{ij} = \mu + \gamma_i + \alpha_j + \varepsilon_{ij}\]

dove \(Y\) è la produzione nel blocco \(i\) e con il diserbo \(j\), \(\mu\) è l’intercetta, \(\gamma\) è l’effetto del blocco \(i\), \(\alpha\) è l’effetto del trattamento \(j\) e \(\varepsilon\) è l’errore sperimentale per ogni singola parcella, che si assume normalmente distribuito, con media 0 e deviazione standard \(\sigma\). Per gli usuali problemi di stimabilità, poniamo i vincoli \(\gamma_1 = 0\) e \(\alpha_1 = 0\) (vincolo sul trattamento), in modo che \(\mu\) rappresenta il valore atteso nel primo blocco e per il primo trattamento in ordine alfabetico. In totale, vi sono 19 parametri da stimare, più \(\sigma\).

9.3 Stima dei parametri

9.3.1 Coefficienti del modello

La stima dei parametri viene eseguita con il metodo dei minimi quadrati. In questo caso l’esperimento è completamente bilanciato e la stima potrebbe essere fatta banalmente, considerando i valori osservati e le medie aritmetiche per gruppo. Infatti, se prendiamo la media generale, la media del primo blocco e la media del primo erbicida riportate nel codice sovrastante, vediamo che il primo erbicida, rispetto alla media generale, ha determinato un aumento di -44.729375 unità. D’altra parte, il primo blocco ha comportanto un aumento di -8.848125 unità. Di conseguenza, il valore atteso nel primo trattamento e nel primo blocco dovrebbe essere pari a:

\[ \bar{Y}_{11} = \mu = 74.56687 - 44.72937 - 8.848125 = 20.98937\]

Analogamente, considerando che la media del secondo blocco è di 2.72625 unità più alta della media generale, il valore atteso per il primo trattamento nel secondo blocco è pari a:

\[ \bar{Y}_{12} = 74.56687 - 44.72937 + 2.72625 = 32.56375\] Di conseguenza:

\[ \gamma_2 = 32.56375 - 20.98937 = 11.57438\]

Ancora, considerando che la media del secondo trattamento è di 14.403125 unità più alta della media generale, il valore atteso per il secondo trattamento nel primo blocco è pari a:

\[ \bar{Y}_{21} = 74.56687 + 14.40313 - 8.848125 = 80.12188\]

Di conseguenza:

\[ \alpha_2 = 80.12188 - 20.98937 = 59.13251\]

Ovviamente, continuare questi calcoli è abbastanza tedioso e quindi eseguiamo il ‘model fitting’ con R, osservando che esso produce gli stessi risultati già calcolati a mano.

options(width = 170)
mod <- lm(Yield ~ factor(Block) + factor(Herbicide), data = rimsulfuron)
summary(mod)
## 
## Call:
## lm(formula = Yield ~ factor(Block) + factor(Herbicide), data = rimsulfuron)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.539  -8.740   1.102   6.209  35.406 
## 
## Coefficients:
##                                                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                                                   20.989      6.886   3.048 0.003847 ** 
## factor(Block)2                                                11.574      4.468   2.590 0.012874 *  
## factor(Block)3                                                17.469      4.468   3.910 0.000309 ***
## factor(Block)4                                                 6.349      4.468   1.421 0.162206    
## factor(Herbicide)Hand-Weeded                                  59.133      8.936   6.617 3.78e-08 ***
## factor(Herbicide)Metolachlor + terbuthylazine (pre)           17.163      8.936   1.921 0.061144 .  
## factor(Herbicide)Pendimethalin (post) + rimsuulfuron (post)   66.795      8.936   7.474 2.03e-09 ***
## factor(Herbicide)Pendimethalin (pre) + rimsulfuron (post)     53.198      8.936   5.953 3.67e-07 ***
## factor(Herbicide)Rimsulfuron (40)                             65.555      8.936   7.336 3.25e-09 ***
## factor(Herbicide)Rimsulfuron (45)                             61.728      8.936   6.907 1.40e-08 ***
## factor(Herbicide)Rimsulfuron (50)                             68.130      8.936   7.624 1.22e-09 ***
## factor(Herbicide)Rimsulfuron (50+30 split)                    59.413      8.936   6.648 3.39e-08 ***
## factor(Herbicide)Rimsulfuron (60)                             58.540      8.936   6.551 4.74e-08 ***
## factor(Herbicide)Rimsulfuron + Atred                          67.148      8.936   7.514 1.77e-09 ***
## factor(Herbicide)Rimsulfuron + hoeing                         64.065      8.936   7.169 5.73e-09 ***
## factor(Herbicide)Rimsulfuron + thyfensulfuron                 54.832      8.936   6.136 1.96e-07 ***
## factor(Herbicide)Thifensulfuron                               22.075      8.936   2.470 0.017359 *  
## factor(Herbicide)Unweeded 1                                   -7.982      8.936  -0.893 0.376473    
## factor(Herbicide)Unweeded 2                                    5.880      8.936   0.658 0.513902    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.64 on 45 degrees of freedom
## Multiple R-squared:  0.8664, Adjusted R-squared:  0.8129 
## F-statistic: 16.21 on 18 and 45 DF,  p-value: 4.916e-14

9.3.2 Residui e devianze

Nel capitolo 7 abbiamo visto come si calcolano i residui, sia a mano che con R. Abbiamo anche visto che la devianza dei residui è la somma dei loro quadrati:

RSS <- sum( residuals(mod)^2 )
RSS
## [1] 7187.348

Da questa possiamo ottenere la deviazione standard (\(\sigma\)), considerando che i gradi di libertà si calcolano partendo dallo stesso punto da cui siamo partiti per l’ANOVA ad una via: abbiamo 16 gruppi con quattro repliche, per cui la devianza, entro ogni gruppo, ha tre gradi di libertà. In totale abbiamo quindi \(16 \times 3 = 48\) gradi di libertà anche se non dobbiamo dimenticare che, le repliche di ogni gruppo non differiscono solo per motivi casuali, ma anche perché appartengono a blocchi diversi. Abbiamo quattro blocchi, quindi tre gradi di libertà, che vanno dedotti dai 48 appena calcolati. Pertanto:

sigma <- sqrt(RSS/45)
sigma
## [1] 12.63799

Più facilmente:

summary(mod)$sigma
## [1] 12.63799

Da \(\sigma\) possiamo ottenere SEM e SED, anche se questo calcolo ve lo lascio per esercizio.

9.4 Scomposizione della varianza

La scomposizione della varianza è sequenziale ed analoga a quella che abbiamo operato per l’ANOVA ad una via; tuttavia, dobbiamo tener presente che, in questo caso, la devianza totale delle osservazione deve essere decomposta in tre quote: una dovuta al trattamento, una dovuta al blocco ed una dovuta agli effetti stocastici.

La devianza totale dei dati, con R, potrebbe essere ottenuta considerando un modello lineare in cui esiste solo l’intercetta, il che equivale a dire che i residui rappresentano gli scostamenti rispetto alla media generale. La devianza dei residui (somma dei quadrati) è quindi la devianza totale delle osservazioni:

mod1 <- lm(Yield ~ 1, data = rimsulfuron)
RSS1 <- sum( residuals(mod1)^2 )
RSS1
## [1] 53779.07

In seconda battuta, inseriamo l’effetto del blocco:

mod2 <- lm(Yield ~ factor(Block), data = rimsulfuron)
RSS2 <- sum( residuals(mod2)^2 )
RSS2
## [1] 51118.58

Vediamo che la devianza del residuo è calata di RSS1 - RSS2 unità, che corrispondono all’effetto del blocco e alla sua introduzione nel modello.

Infine, inseriamo anche l’effetto del trattamento, tornando così al modello completo che abbiamo utilizzato più sopra. Considerando i paragrafi precedenti, notiamo che l’inserimento dell’effetto del trattamento ha determinato un calo della devianza dei residui pari RSS2 - RSS, che corrisponde appunto all’effetto del diserbo. Insomma, ogni effetto introdotto nel modello ha determinato un decremento della variabilità non spiegata, quantitativamente pari alla variabilità attribuibile all’effetto stesso. Alla fine del processo rimane comunque un certo residuo (RSS) non spiegato, corrispondente a \(\varepsilon\).

Ci chiediamo se gli effetti attribuibili al blocco e al trattamento sono significativamente più grandi del residuo. Sappiamo già di non poter confrontare le devianze, ma possiamo calcolare e confrontare con un test di F le relative varianze. Basta tener conto che i gradi di libertà dei blocchi e dei trattamenti sono rispettivamente 3 e 15, cioè il numero dei blocchi meno uno e il numero dei trattamenti meno uno.

La tabella ANOVA può essere facilmente ottenuta con R:

anova(mod)
## Analysis of Variance Table
## 
## Response: Yield
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## factor(Block)      3   2660  886.83  5.5524  0.002496 ** 
## factor(Herbicide) 15  43931 2928.75 18.3369 2.329e-14 ***
## Residuals         45   7187  159.72                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ovviamente, prima di considerare questa tabella dovremo preoccuparci del fatto che le assunzioni di base siano rispettate, cosa che possiamo facilmente verificare con un’analisi grafica dei residui, utilizzando il codice sottostante. L’output è visibile in Figura 9.1.

par(mfrow=c(1,2))
plot(mod, which = 1)
plot(mod, which = 2)
Analisi grafica dei residui per la prova di confronto erbicida

Figure 9.1: Analisi grafica dei residui per la prova di confronto erbicida

Dopo esserci rassicurati su questo importante aspetto, possiamo vedere che abbiamo due ipotesi nulle da testare (effetto del trattamento non significativo ed effetto del blocco non significativo), che possono essere entrambe rifiutate per P < 0.05.

Da questo punto in avanti, l’analisi procede come usuale, calcolando le medie marginali attese ed, eventualmente, confrontandole tra loro con una procedura di confronto multiplo, come descritto nei capitoli precedenti. Tener presente che, in questo esperimento, abbiamo 16 trattamenti, cioè \(16 \times 15 / 2 = 120\) confronti; di conseguenza, può essere opportuno operare la correzione per la molteplicità. Inoltre, dato che il trattamento più interessante è quello che rende massima la produzione, sarà opportuno ordinare le medie in senso decrescente, utilizzando l’argomento ‘reverse = T’.

library(emmeans)
medie <- emmeans(mod, ~factor(Herbicide))
multcomp::cld(medie, Letters = LETTERS, reverse = T)
##  Herbicide                                  emmean   SE df lower.CL upper.CL .group
##  Rimsulfuron (50)                             98.0 6.32 45    85.24    110.7  A    
##  Rimsulfuron + Atred                          97.0 6.32 45    84.26    109.7  A    
##  Pendimethalin (post) + rimsuulfuron (post)   96.6 6.32 45    83.91    109.4  A    
##  Rimsulfuron (40)                             95.4 6.32 45    82.67    108.1  A    
##  Rimsulfuron + hoeing                         93.9 6.32 45    81.18    106.6  A    
##  Rimsulfuron (45)                             91.6 6.32 45    78.84    104.3  A    
##  Rimsulfuron (50+30 split)                    89.2 6.32 45    76.52    102.0  A    
##  Hand-Weeded                                  89.0 6.32 45    76.24    101.7  A    
##  Rimsulfuron (60)                             88.4 6.32 45    75.65    101.1  A    
##  Rimsulfuron + thyfensulfuron                 84.7 6.32 45    71.94     97.4  A    
##  Pendimethalin (pre) + rimsulfuron (post)     83.0 6.32 45    70.31     95.8  AB   
##  Thifensulfuron                               51.9 6.32 45    39.19     64.6   BC  
##  Metolachlor + terbuthylazine (pre)           47.0 6.32 45    34.27     59.7    C  
##  Unweeded 2                                   35.7 6.32 45    22.99     48.4    C  
##  Alachlor + terbuthylazine                    29.8 6.32 45    17.11     42.6    C  
##  Unweeded 1                                   21.9 6.32 45     9.13     34.6    C  
## 
## Results are averaged over the levels of: Block 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 16 estimates 
## significance level used: alpha = 0.05

Vi lascio il commento dei risultati come esercizio.

9.5 Disegni a quadrato latino

A volte i fattori di blocco sono più di uno e danno origine ad un disegno sperimentale detto ‘quadrato latino’ di cui abbiamo parlato nel capitolo 3. Qui, forniremo un esempio tratto dalla pratica sperimentale industriale.

9.6 Caso studio: confronto tra metodi costruttivi

Immaginiamo di voler studiare il tempo necessario per costruire un componente elettronico, utilizzando quattro metodi diversi. E’ evidente che il tempo di costruzione sarà influenzato dalla perizia del tecnico e, per questo, utilizziamo quattro tecnici diversi, ad ognuno dei quali facciamo utilizzare tutti e quattro i metodi. Un esperimento così disegnato sarebbe a blocchi randomizzati, con il tecnico che fa da blocco per i trattamenti. Tuttavia, dobbiamo anche riconoscere che i quattro tecnici saranno via via meno efficienti, e quindi il metodo che utilizzeranno per primo sarà avvantaggiato, mentre quello che utilizzeranno per ultimo sarà svantaggiato. E’vero che i metodi sono assegnati in ordine random ad ogni tecnico, ma non si può comunque evitare che un metodo venga avvantaggiato rispetto ad un altro, perché, ad esempio, non viene mai ad occupare l’ultima posizione (o meglio, l’ultimo turno).

Per evitare questo problema imponiamo un vincolo ulteriore alla randizzazione e facciamo in modo che ogni metodo occupi tutte e quattro i turni, in tecnici diversi. Il disegno è quindi a quadrato latino.

Il dataset dei risultati è disponibile su gitHub:

path1 <- "https://raw.githubusercontent.com/OnofriAndreaPG/"
path2 <- "aomisc/master/data/"
fileName <- "Technicians.csv"
filePath <- paste(path1, path2, fileName, sep = "")
dataset <- read.csv(filePath, header=T)
dataset
##    Shift Technician Method Time
## 1      I     Andrew      C   90
## 2     II     Andrew      B   90
## 3    III     Andrew      A   89
## 4     IV     Andrew      D  104
## 5      I       Anna      D   96
## 6     II       Anna      C   91
## 7    III       Anna      B   97
## 8     IV       Anna      A  100
## 9      I    Michael      A   84
## 10    II    Michael      D   96
## 11   III    Michael      C   98
## 12    IV    Michael      B  104
## 13     I      Sarah      B   88
## 14    II      Sarah      A   88
## 15   III      Sarah      D   98
## 16    IV      Sarah      C  106

9.7 Definizione di un modello lineare

In questo caso abbiamo un trattamento (metodo) e due effetti ‘blocco’ (tecnico e turno) da includere nel modello, che può essere così definito:

\[Y_{ijk} = \mu + \gamma_k + \beta_j + \alpha_i + \varepsilon_{ijk}\]

dove \(\mu\) è l’intercetta, \(\gamma\) è l’effetto del turno k, \(\beta\) è l’effetto del tecnico j e \(\alpha\) è l’effetto del metodo i. L’elemento \(\varepsilon_{ijk}\) rappresenta la componente random individuale, di ogni osservazione e si assume normalmente distribuita, con media 0 e deviazione standard \(\sigma\).

Avendo già illustrato il processo di stima dei parametri e di scomposizione della varianza e quindi utilizziamo subito R per il ‘model fitting’:

mod <- lm(Time ~ Method + Technician
          + Shift, data = dataset)

Verifichiamo il rispetto delle assunzioni di base, con l’analisi grafica dei residui, riportata in Figura 9.2 (il codice è analogo a quello fornito più sopra).

Analisi grafica dei residui per la prova di confronto tra metodi costruttivi

Figure 9.2: Analisi grafica dei residui per la prova di confronto tra metodi costruttivi

Non essendovi evidenti problemi, valutiamo la significatività degli effetti nel modello, analogamente a quanto abbiamo fatto nel caso dell’ANOVA a blocchi randomizzati. L’unica differenza sta nel fatto che, nei disegni a quadrato latino, vi sono tre effetti da testare, anche se l’unico ad avere una certa rilevanza è l’effetto del metodo di lavoro.

anova(mod)
## Analysis of Variance Table
## 
## Response: Time
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Method      3 145.69  48.563 12.7377 0.0051808 ** 
## Technician  3  17.19   5.729  1.5027 0.3065491    
## Shift       3 467.19 155.729 40.8470 0.0002185 ***
## Residuals   6  22.87   3.812                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vediamo che esiste una differenza significativa tra i metodi e l’ipotesi nulla può essere rifiutata con una bassissima probabilità di errore di prima specie.

Ovviamente, dopo aver eseguito un’ANOVA a blocchi randomizzati o a quadrato latino, andremo eventualmente ad eseguire un test di confronto multiplo, seguendo le raccomandazioni esposte nel capitolo precedente. Anche questa parte ve la lascio per esercizio.