Capitolo 6 Modelli ANOVA ad una via

Nel capitolo 4 abbiamo già parlato di come assumiamo che risultati di un esperimento si generino attraverso un doppio meccanismo deterministico e stocastico, rappresentabile attraverso un modello matematico contenento uno o più elementi casuali, descritti da una funzione di densità, solitamente gaussiana. Abbiamo già visto che i modelli più utilizzati sono lineari e, tra questi, il gruppo più diffuso, almeno in ambito biologico-agrario, è quello dei modelli ANOVA.

La nomenclatura è impropria; infatti, con il termine ANOVA, si intende un’operazione di scomposizione della varianza, inventata da Fisher negli anni ‘30 del ’900, che costituisce ancora oggi una delle più fondamentali tecniche di analisi dei dati sperimentali. Tuttavia, ci sono alcuni modelli che, più naturalmente di altri, sono connessi all’esecuzione dell’ANOVA fisheriana. Si tratta di modelli lineari nei quali la/le variabile/i indipendenti, che descrivono gli stimoli sperimentali, sono nominali (categoriche) e, nella letteratura anglosassone, prendono il nome di ’factors’ (fattori sperimentali).

Questi modelli rappresentano il punto di ingresso nell’analisi dei dati e la gran parte della letteratura scientifica è basata proprio su questi modelli. Pertanto, ad essi dedicheremo ampio spazio nel seguito di questo libro.

6.1 Caso-studio: confronto tra erbicidi in vaso

Abbiamo eseguito un esperimento in vaso, nel quale abbiamo utilizzato quattro trattamenti erbicidi:

  1. Metribuzin
  2. Rimsulfuron
  3. Metribuzin + rimsulfuron
  4. Testimone non trattato

Lo scopo dell’esperimento era quello di verificare se la miscela metribuzin e rimsulfuron è più efficace dei due componenti utilizzati separatamente. L’esperimento era disegnato a randomizzazione completa con quattro repliche e prevedeva il rilievo della biomassa presente su ogni vaso, tre settimane dopo il trattamento: un più basso valore di biomassa implica un miglior effetto del trattamento.

I risultati di questo esperimento sono riportati nel dataset ‘mixture.csv’, che è disponibile su gitHub. Per prima cosa, carichiamo il file.

str1 <- "https://raw.githubusercontent.com/OnofriAndreaPG"
str2 <- "/agroBioData/master/mixture.csv"
pathData <- paste(str1, str2, sep = "")

dataset <- read.csv(pathData, header = T)
head(dataset)
##             Treat Weight
## 1        Unweeded  24.62
## 2        Unweeded  30.94
## 3        Unweeded  24.02
## 4        Unweeded  27.51
## 5 Metribuzin__348  15.20
## 6 Metribuzin__348   4.38

Il dataset è organizzato come un database, nel quale ogni riga contraddistingue un’unità sperimentale (record) e ogni colonna rappresenta una caratteristica del record (campo). In questo caso abbiamo due colonne: una per la variabile indipendente (fattore sperimentale ‘Treat’), che ci dice con quale erbicida è stata trattata ogni unità sperimentale, ed una per la variabile dipendente (‘Weight’).

6.2 Descrizione del dataset

La prima analisi dei dati consiste nella valutazione descrittiva del dataset. In particolare, calcoliamo:

  1. le medie per ogni tesi sperimentale
  2. Le deviazioni standard per ogni tesi sperimentale

Utilizziamo la funzione ‘tapply’ in R.

medie <- tapply(dataset$Weight, dataset$Treat, mean)
SDs <- tapply(dataset$Weight, dataset$Treat, sd)
descrit <- cbind(medie, SDs)
descrit
##                   medie      SDs
## Metribuzin__348  9.1750 4.699089
## Mixture_378      5.1275 2.288557
## Rimsulfuron_30  16.8600 4.902353
## Unweeded        26.7725 3.168673

Che cosa ci dice questa tabella, in base agli obiettivi dell’esperimento?

Ci suggerisce che:

  1. La miscela sembra leggermente più efficace dei prodotti singoli
  2. Esiste una certa variabilità (errore sperimentale), che impedisce un giudizio certo
  3. La variabilità è abbastanza simile per tutti i trattamenti

6.3 Definizione di un modello lineare

Sappiamo che, dietro ai campioni osservati, vi sono una o più popolazioni di riferimento, alle quali è rivolto il nostro interesse: abbiamo infatti osservato \(m_1\), \(m_2\), \(m_3\) ed \(m_4\) (le medie dei campioni), ma siamo effettivamente interessati a conoscere \(\mu_1\), \(\mu_2\), \(\mu_3\) e \(\mu_4\), cioè le medie delle popolazioni da cui i campioni sono tratti. E siamo anche interessati a capire se la differenza è sufficientemente ampia da poter ritenere che non sia di natura puramente casuale.

Per descrivere le osservazioni possiamo utilizzare un modello lineare del tipo:

\[Y_i = \mu + \alpha_j + \varepsilon_i\]

Questo modello impone che le osservazioni \(Y\) derivino da una valore \(\mu\), detto intercetta, a cui si aggiunge l’effetto del trattamento \(\alpha_j\) e l’effetto stocastico \(\varepsilon\), che è un elemento che influenza ogni singola osservazione \(i\).

Il valore atteso per un soggetto, in assenza di errore sperimentale, è:

\[\bar{Y_i} = \mu + \alpha_j\]

È abbastanza facile intuire che, se non ci fosse l’errore sperimentale, un soggetto dovrebbe avere un valore pari alla media del gruppo di cui fa parte, per cui:

\[\bar{Y_i} = \mu_j = \mu + \alpha_j\]

6.4 Parametrizzazione del modello

Gli elementi \(\mu\) ed \(\alpha\) hanno un significato diverso a seconda della parametrizzazione del modello. Infatti, dobbiamo tener presente che, con il modello lineare sopra riportato, ci stiamo proponendo di ottenere la quantità \(Y_i\) come somma di tre valori ed, in questi termini, il problema è indeterminato. Infatti esistono infinite triplette di valori che, sommati, possono restituire \(Y_i\), qualunque esso sia. Per questo motivo dobbiamo porre un qualche vincolo su uno degli elementi in gioco. Esistono diversi tipi di vincoli, ma, in questo testo, ne tratteremo solo due, quelli più rilevanti.

6.4.1 Vincolo sul trattamento

Un vincolo molto usato è \(\alpha_1 = 0\). Di conseguenza, risulta che:

\[ \left\{ {\begin{array}{l} \mu_1 = \mu + \alpha_1 = \mu + 0\\ \mu_2 = \mu + \alpha_2 \\ \mu_3 = \mu + \alpha_3 \\ \mu_4 = \mu + \alpha_4 \end{array}} \right.\]

quindi \(\mu\) è la media del primo trattamento, usualmente inteso in ordine alfabetico (ma il riferimento può cambiare, a seconda delle esigenze). Gli altri valori \(\alpha_j\), sono differenze tra le medie dei gruppi da 2 a 4 con la media del gruppo 1.

In generale, adottando il vincolo sul trattamento, i parametri sono medie, valori attesi oppure differenze tra medie e/o differenze tra valori attesi.

6.4.2 Vincolo sulla somma

Un altro possibile vincolo consiste nell’imporre che la somma degli \(\alpha_j\) sia pari a 0. In questo caso, abbiamo quattro gruppi, quindi quattro medie. Se prendiamo l’espressione precedente e sommiamo membro a membro, otteniamo:

\[ \mu_1 + \mu_2 + \mu_3 + \mu_4 = 4 \mu + \sum{\alpha_j}\]

Quindi, con il vincolo sulla somma:

\[ \mu_1 + \mu_2 + \mu_3 + \mu_4 = 4 \mu\] Quindi:

\[\mu = \frac{\mu_1 + \mu_2 + \mu_3 + \mu_4}{4} \]

cioè \(\mu\) è la media generale e \(\alpha_j\) rappresentano gli scostamenti di ogni trattamento rispetto alla media generale, usualmente definiti effetti dei trattamenti. Se un prodotto è efficace, abbasserà di più il peso delle infestanti e quindi avrà un elevato effetto negativo.

La seconda parametrizzazione è forse più ‘comprensibile’ e, se il disegno è bilanciato, è abbastanza facile stimare i parametri ‘a mano’, utilizzando il metodo dei momenti, a partire dalle medie aritmetiche dei trattamenti.

I software statistici, invece, non utilizzano il metodo dei momenti, ma scelgono i valori dei parametri ai quali corrisponde il minimo valore della somma dei quadrati dei residui (metodo dei minimi quadrati). Solitamente utilizzano di default la parametrizzazione con vincolo sul trattamento, anche se questa impostazione, se necessario, può essere cambiata.

Indipendentemente dalla parametrizzazione prescelta, i valori attesi, e quindi i residui, sono gli stessi; l’errore sperimentale viene assunto gaussiano e omoscedastico, cioè la varianza è unica ed uguale per tutte le tesi:

\[\varepsilon_i \sim N(0, \sigma)\]

6.5 Assunzioni di base

In questa costruzione algebrica abbiamo implicitamente posto alcuni ‘punti fermi’, detti assunzioni di base, che sono i seguenti:

  1. la componente deterministica è lineare e additiva
  2. non vi sono altri effetti, se non il trattamento e l’errore, che è puramente stocastico, senza componenti sistematiche
  3. gli errori sono campionati in modo indipendente, da una distribuzione normale, con media 0 e deviazione standard \(\sigma\)
  4. le varianze sono omogenee (unico valore di \(\sigma\), comune per tutti i gruppi)

È evidente che il nostro dataset deve conformarsi a queste nostre aspettative, altrimenti il modello è invalido. Ci occuperemo di questa verifica nel prossimo capitolo.

6.6 Stima dei parametri

6.6.1 Coefficienti del modello

Iniziamo il nostro lavoro di inferenza con la stima dei parametri. Dato che il disegno è bilanciato (stesso numero di repliche per trattamento), la stima dei parametri può essere fatta a mano, a partire dalle medie aritmetiche dei trattamenti e dalla media generale (metodo dei momenti). Avendo scelto di imporre un vincolo sul trattamento (\(\alpha_1 = 0\)), i valori dei parametri possono essere così ottenuti:

\[ \left\{ {\begin{array}{l} \mu = 9.175\\ \alpha_2 = \mu_2 - \mu = 5.1275 - 9.1750 = - 4.0475\\ \alpha_3 = \mu_3 - \mu = 16.86 - 9.1750 = 7.685\\ \alpha_4 = \mu_4 - \mu = 26.7725 - 9.1750 = 17.5975\\ \end{array}} \right.\]

6.6.2 Residui

Abbiamo già visto che i valori attesi (\(\bar{Y_i}\)) sono pari, per ogni osservazione, alla somma tra \(\mu\) e il valor \(\alpha\) relativo al gruppo di cui l’osservazione fa parte. In questo caso, il valore atteso coincide con la media di ogni gruppo e, pertanto, i residui (\(\varepsilon_i = Y - \bar{Y_i}\)) possono essere calcolati facilmente, come indicato in Tabella 6.1.

Table 6.1: Tabella dei dati osservati, dei valori attesi, dei residui e dei loro quadrati
Erbicida Weight Attesi Residui Residui^2
Unweeded 24.62 26.773 -2.153 4.633
Unweeded 30.94 26.773 4.167 17.368
Unweeded 24.02 26.772 -2.752 7.576
Unweeded 27.51 26.772 0.738 0.544
Metribuzin__348 15.20 9.175 6.025 36.301
Metribuzin__348 4.38 9.175 -4.795 22.992
Metribuzin__348 10.32 9.175 1.145 1.311
Metribuzin__348 6.80 9.175 -2.375 5.641
Mixture_378 6.14 5.127 1.013 1.025
Mixture_378 1.95 5.127 -3.177 10.097
Mixture_378 7.27 5.127 2.143 4.590
Mixture_378 5.15 5.127 0.023 0.001
Rimsulfuron_30 10.50 16.860 -6.360 40.450
Rimsulfuron_30 20.70 16.860 3.840 14.746
Rimsulfuron_30 20.74 16.860 3.880 15.054
Rimsulfuron_30 15.50 16.860 -1.360 1.850

6.6.3 Stima di \(\sigma\)

Abbiamo detto che \(\sigma\) è la deviazione standard dei residui, che tuttavia deve essere calcolata facendo attenzione al numero dei gradi di libertà. Partiamo dalla devianza dei residui, che possiamo ottenere come somma dei quadrati riportati nell’ultima colonna a destra in Tabella 6.1. Otteniamo il valore RSS = 184.17745.

I residui costituiscono lo scostamento di ogni dato rispetto alla media del trattamento e, di conseguenza, per ogni gruppo la loro somma deve essere 0. Quindi per ogni gruppo vi sono 3 gradi di libertà, cioè 12 gradi di libertà in totale (3 \(\times\) 4, cioè il numero dei trattamenti per il numero delle repliche meno una). Ne consegue che la varianza dei residui è pari a:

\[MS_{e} = \frac{184.178}{12} = 15.348\]

La deviazione standard \(\sigma\) è:

\[ \sqrt{15.348} = 3.9177\]

6.6.4 SEM e SED

Abbiamo visto che la varianza d’errore è pari a 15.348 e pertiene ad ogni singola osservazione effettuata durante l’esperimento. Questa osservazione ci può aiutare a costruire una banda di inferenza attorno alle medie stimate; infatti, noi abbiamo osservato \(m_1\), \(m_2\), \(m_3\) ed \(m_4\) e con la nostra stima puntuale, abbiamo assunto che i valori osservati coincidessero rispettivamente con \(\mu_1\), \(\mu_2\), \(\mu_3\) e \(\mu_4\). Analogamente a quanto abbiamo visto in un capitolo precedente, possiamo calcolare l’incertezza associata alle nostre stime attraverso l’Errore Standard di una Media (SEM), pari a:

\[SEM = \sqrt{ \frac{MS_e}{n} } = \frac{3.9177}{\sqrt{4}}\]

Oltre che di una media, spesso siamo interessati anche a conoscere la varianza della differenza di due medie. Dato che la differenza di variabili casuali ha una varianza pari alla somma delle varianza delle due variabili originali, possiamo scrivere che:

\[SED = \sqrt{ MS_{media1} + MS_{media2} } = \sqrt{ 2 \cdot \frac{MS_e}{n} } = \sqrt{2} \cdot \frac{3.9177}{\sqrt{4}} = \sqrt{2} \cdot SEM\]

6.7 Scomposizione della varianza

La scomposizione della varianza, il cui significato sarà più chiaro in seguito, è la vera e propria ANOVA fisheriana. Può essere eseguita seguendo un metodo sequenziale, che è molto semplice e perfettamente valido quando i dati sono bilanciati, come in questo caso.

Il primo elemento da stimare è la devianza totale dei dati, cioè la somma dei quadrati degli scarti di ogni dato rispetto alla media generale:

\[\begin{array}{c} SS = \left(24.62 - 14.48375\right)^2 + \left(30.94 - 14.48375\right)^2 + ... \\ ... + \left(15.50 - 14.48375\right)^2 = 1273.706 \end{array}\]

Questa devianza esprime la variabilità totale tra un dato e l’altro, sia quella dovuta al trattamento, sia quella puramente stocastica.

Abbiamo già calcolato la devianza dei residui (RSS = 184.17745), che quantifica la variabilità dei dati dovuta ad effetti puramente casuali, ma non al trattamento sperimentale. Infatti, la variabilità di un dato rispetto alla media del gruppo a cui appartiene non può essere dovuta al trattamento, in quanto i soggetti dello stesso gruppo sono trattati allo stesso modo.

A questo punto ci chiediamo: qual è la devianza prodotta dal trattamento sperimentale? La devianza totale rappresenta la variabilità totale dei dati (trattamento + effetti stocastici), la devianza residua rappresenta la sola variabilità stocastica, di conseguenza, la differenza di queste due quantità rappresenta la devianza dovuta al trattamento sperimentale:

\[SS_t = SS - RSS = 1273.706 + 184.1774 = 1089.529\]

In realtà, oltreché per differenza, la devianza del trattamento potrebbe anche essere ottenuta direttamente. Dobbiamo tener conto che la media generale delle osservazioni è pari a 14.48375; se le medie dei gruppi differiscono dalla media generale, ciò potrebbe essere dovuto all’effetto del trattamento, che è diverso in ogni gruppo. Quindi, la devianza del trattamento può essere calcolata come il quadruplo (abbiamo quattro repliche) della devianza tra le medie dei trattamenti:

\[{\begin{array}{l} SS_t = 4 \times \left[ \left(9.1750 - 14.48375\right)^2 + \left(5.1275 - 14.48375\right)^2 + \right. \\ + \left. \left(16.86 - 14.48375\right)^2 + \left(26.7725 - 14.48375\right)^2 \right] = 1089.529 \end{array}} \]

A questo punto, ricordiamo che il nostro obiettivo è stabilire se il trattamento ha avuto un effetto significativo, cioè se l’effetto da esso prodotto si distingue dal ‘rumore di fondo’, ossia dall’elemento casuale rappresentato dal residuo. Tuttavia, le devianza del trattamento e del residuo non possono essere confrontate direttamente, in quanto hanno un numero diverso di gradi di libertà.

In particolare, la devianza del trattamento ha tre gradi di libertà (tre scarti liberi; in genere numero dei trattamenti - 1) e quindi la relativa varianza è:

\[MS_t = \frac{1089.529}{3} = 363.1762\]

La varianza del residuo è invece pari a 15.348 e, come abbiamo visto più sopra, ha 12 gradi di libertà. Le due varianze così ottenute possono essere confrontate, all’interno di una procedura di test d’ipotesi, che vedremo di seguito.

Prima, però, vale la pena di sottolineare come la procedura precedente ci abbia permesso di suddividere la variabilità totale delle osservazioni in due componenti, una dovuta al trattamento e una dovuta all’errore sperimentale. Per questo motivo, parliamo di scomposizione della varianza (variance partitioning) o analisi della varianza (ANalysis Of VAriance = ANOVA), che è lo strumento più utilizzato nella sperimentazione biologica e che dobbiamo totalmente all’inventiva di Ronald Fisher.

6.8 Test d’ipotesi

Abbiamo descritto il nostro esperimento e ne abbiamo individuato le caratteristiche rilevanti, stimando i parametri che meglio le descrivono (effetti dei trattamenti e varianza). Come già detto, dobbiamo chiederci se i dati rispettano le assunzioni di base del modello, ma di questo parleremo in una lezione a parte.

Ora, il nostro scopo è capire se il trattamento sperimentale abbia prodotto un effetto rilevante, maggiore di quello prodotto da altri elementi puramente casuali (‘rumore di fondo’).

L’ipotesi nulla è che il trattamento non abbia avuto effetto, cioè che:

\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\]

In altri termini, se l’ipotesi nulla è vera, i quattro campioni sono estratti da quattro popolazioni identiche, o meglio, dalla stessa popolazione. Ciò può essere anche declinato come:

\[H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 0\]

Una statistica rilevante per testare questa ipotesi è data dal rapporto tra la varianza del trattamento e quella dell’errore:

\[F = \frac{MS_t}{MS_e} = \frac{363.18}{15.348} = 23.663\]

E’ evidente che se il trattamento non fosse stato efficace, non dovrebbe aver prodotto una variabilità di molto superiore a quella dell’errore (quindi F = 1). In questo caso la variabilità prodotta dal trattamento è stata oltre 23 volte superiore a quella dell’errore. Delle due l’una: o il trattamento è stato efficace oppure io sono stato particolarmente sfortunato e, nell’organizzare questo esperimento, si è verificato un caso particolarmente raro.

Ci chiediamo: “se l’ipotesi nulla è vera, qual è la probabilità di osservare un valore di F così alto o più alto?”. In altre parole, “qual è la sampling distribution per F?”. Potremmo costruire questa distribuzione empiricamente, attraverso una simulazione Monte Carlo.

Se assumiamo che l’ipotesi nulla è vera, allora i dati dovrebbero essere campionati da una distribuzione gaussiana con media pari a 14.48375 e deviazione standard pari a 3.9177 (vedi sopra). Possiamo quindi utilizzare un generatore di numeri casuali gaussiani, per ottenere un dataset, per il quale l’ipotesi nulla è certamente vera e sottoporlo alla scomposizine della varianza, per calcolare un valore di F.

Se ripetiamo questo processo, ad esempio, 100’000 volte, otteniamo 100’000 dataset, nei quali non vi è effetto del trattamento e altrettanti valori F in situazioni nelle quali l’ipotesi nulla è vera. Otteniamo quindi una sampling distribution empirica per F e possiamo valutare quanti valori sono superiori a 23.66.

Per chi voglia eseguire questa simulazione, il codice da utilizzare è riportato nell’appendice finale di questo capitolo. Diciamo solo che il valore minimo di F è stato 0.00019, il massimo 32.87, la media è stata 1.19 e la mediana 0.84. Tra tutti i 100’000 valori, ne abbiamo trovati solo due pari o superiori a quello osservato, il che vuol dire che la probabilità che l’ipotesi nulla sia vera con F = 23.66 è pari a 2 \(\times 10^{-5}\).

La sampling distribution (opportunamente discretizzata) è riportata in figura 6.1. Si tratta di una distribuzione chiaramente non normale, ma assimilabile alla distribuzione F di Fisher (in realtà l’invenzione è di Snedecor, allievo di Fisher), con 3 gradi di libertà al numeratore e 12 al denominatore (in blue in figura).

Sampling distribution empirica per F, con l'ipotesi nulla vera, per l'esempio

Figure 6.1: Sampling distribution empirica per F, con l’ipotesi nulla vera, per l’esempio

Di conseguenza, possiamo evitare la simulazione Monte Carlo ed utilizzare la F di Fisher per calcolare la probabilità di ottenere un valore di F altrettanto estremo o più estremo del nostro. Ad esempio, in R, possiamo utilizzare la funzione:

pf(23.663, 3, 12, lower.tail = F)
## [1] 2.508789e-05

che porta ad un risultato molto simile a quello ottenuto con la simulazione di Monte Carlo. Insomma, in assenza di un effetto del trattamento (quindi per il solo effetto del caso), se ripetiamo l’esperimento infinite volte, abbiamo una probabilità molto bassa che si produca un valore di F altrettanto alto o più alto di quello da noi osservato.

Di conseguenza, se rifiutiamo l’ipotesi nulla di assenza di effetto del trattamento e accettiamo l’ipotesi alternativa (il trattamento ha avuto effetto significativo) ci portiamo dietro un rischio di errore estremamente piccolo, comunque molto al disotto della soglia prefissata del 5%.

6.9 Operazioni con R

La stima dei parametri di un modello lineare, in R, può essere eseguita in modo molto banale, utilizzando la funzione ‘lm()’ Il codice è il seguente:

mod <- lm(Weight ~ factor(Treat), data = dataset)

Il primo argomento rappresenta il modello lineare che vogliamo adattare ai dati: l’incusione dell’intercetta, codificata con il carattere ‘1’ è opzionale (‘Weight ~ 1 + Treat’ e ‘Weight ~ Treat’ sono equivalenti), mentre il termine stocastico \(\epsilon\) viene incluso di default e non deve essere indicato. Il termine ‘factor’ sta a significare che la variabile ‘Treat’ è un fattore sperimentale; questa specifica è opzionale quando la variabile è di tipo ‘carattere’ (come in questo caso), ma è fondamentale quando abbiamo a che fare con una variabile a codifica numerica.

Una volta operato l’adattamento, l’oggetto ‘mod’ contiene tutti i risultati, che possiamo ottenere con opportuni ‘estrattori’. Ad esempio, la funzione ‘summary()’ restituisce le stime dei parametri:

summary(mod)
## 
## Call:
## lm(formula = Weight ~ factor(Treat), data = dataset)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.360 -2.469  0.380  2.567  6.025 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.175      1.959   4.684 0.000529 ***
## factor(Treat)Mixture_378      -4.047      2.770  -1.461 0.169679    
## factor(Treat)Rimsulfuron_30    7.685      2.770   2.774 0.016832 *  
## factor(Treat)Unweeded         17.598      2.770   6.352 3.65e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.918 on 12 degrees of freedom
## Multiple R-squared:  0.8554, Adjusted R-squared:  0.8193 
## F-statistic: 23.66 on 3 and 12 DF,  p-value: 2.509e-05

Notiamo immediatamente che viene utilizzata la parametrizzazione con vincolo sul trattamento, visto che l’intercetta coincide con la media del primo trattamento. I valori di \(\alpha\) non sono altro che differenze tra questa media e tutte le altre medie.

I valori attesi e i residui possono essere ottenuti, in R, applicando due metodi ‘fitted()’ e ‘residuals()’ all’oggetto ottenuto dal fitting:

attesi <- fitted(mod)
epsilon <- residuals(mod)

La devianza del residuo (somma dei quadrati degli scarti) è:

deviance(mod)
## [1] 184.1774

ementre la deviazione standard del residuo (\(\sigma\)), già visualizzabile nell’output della funzione ‘summary()’, può essere estratta come:

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

La tabella finale dell’ANOVA può essere ottenuta in R utilizzando la seguente funzione:

anova(mod)
## Analysis of Variance Table
## 
## Response: Weight
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## factor(Treat)  3 1089.53  363.18  23.663 2.509e-05 ***
## Residuals     12  184.18   15.35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.10 Medie marginali attese

Abbiamo visto che le somme \(\mu + \alpha_j\) restituiscono le medie dei trattamenti e prendono il nome di medie marginali attese. Quando l’esperimento è bilanciato (ugual numero di repliche per ogni trattamento), esse sono uguali alle medie aritmetiche di ogni gruppo, mentre, quando l’esperimento è sbilanciato esso sono diverse, ma costituiscono una stima migliore delle medie delle popolazioni da cui i gruppi sono estratti rispetto alle medie aritmetiche. In R, il metodo più semplice per calcolare le medie marginali attese è quello di utilizzare la funzione ‘emmeans()’, nel package ‘emmeans’ (Lenth, 2016), che deve essere preventivamente installato e caricato in memoria.

#install.packages("emmeans") # Installare il package
library(emmeans) #Caricare il package in memoria
medie <- emmeans(mod, ~Treat)
medie
##  Treat           emmean   SE df lower.CL upper.CL
##  Metribuzin__348   9.18 1.96 12     4.91     13.4
##  Mixture_378       5.13 1.96 12     0.86      9.4
##  Rimsulfuron_30   16.86 1.96 12    12.59     21.1
##  Unweeded         26.77 1.96 12    22.50     31.0
## 
## Confidence level used: 0.95

6.11 Per concludere …

  1. Con l’ANOVA la variabilità totale dei dati viene decomposta in due quote, una attribuibile al trattamento sperimentale ed una inspiegabile (residuo)
  2. L’effetto del trattamento è significativo, se la variabilità che esso provoca è superiore a quella inspiegabile
  3. Il confronto tra varianze viene impostato sotto forma di rapporto (F).
  4. L’ipotesi nulla è che il trattamento NON HA AVUTO effetto significativo
  5. Sotto questa ipotesi nulla, l’F osservato ha una distribuzione F di Fisher
  6. L’ipotesi nulla è rifiutata quando P \(\leq\) 0.05 (livello di protezione arbitrario, ma universalmente accettato)

Ovviamente, è necessario ricordare che tutte le considerazioni fatte finora sono valide se il dataset è conforme alle assunzioni di base per l’ANOVA, per cui bisogna sempre eseguire i necessari controlli, di cui parleremo nel prossimo capitolo.


6.12 Per approfondire un po’…

6.12.1 Simulazione della sampling distribution per F

Di seguito forniamo il codice per simulare la sampling distribution per F, quando l’ipotesi nulla è vera. Il codice impiega un generatore di numeri casuali gaussiano per ottenere dataset nei quali l’ipotesi nulla è vera, perché non vi è effetto del trattamento.

Fvals <- c()
set.seed(1234)
for(i in 1:100000){
  Ysim <- rnorm(16, 14.48375, 3.9177)
  Fvals[i] <- anova ( lm(Ysim ~ factor(Treat), data = dataset))$F[1]
}
min(Fvals)
max(Fvals)
length(Fvals[Fvals > 23.663])

I risultati ottenuti sono quelli esposti più sopra

6.12.2 Altre letture

  1. Faraway, J.J., 2002. Practical regression and Anova using R. http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf, R.
  2. Fisher, Ronald (1918). “Studies in Crop Variation. I. An examination of the yield of dressed grain from Broadbalk” (PDF). Journal of Agricultural Science. 11 (2): 107–135.
  3. Kuehl, R. O., 2000. Design of experiments: statistical principles of research design and analysis. Duxbury Press (CHAPTER 2)
  4. Lenth, R.V., 2016. Least-Squares Means: The R Package lsmeans. Journal of Statistical Software 69. https://doi.org/10.18637/jss.v069.i01