sabato 2 gennaio 2010

Latin squares design in R

Il disegno a quadrati latini (latin square design) permette di analizzare 3 fattori contemporaneamente, così come si farebbe con una ANOVA a 3 fattori.
Questo metodo è stato creato inizialmente in ambito agrario, per avere un metodo di analisi della produttività di un terreno, coltivato in maniera differente. Pertanto userò un esempio prettamente agrario, ma questa analisi può essere applicata anche in altri ambiti.
In un quadrato latino si crea una tabella NxN, con le righe e le colonne rappresentati rispettivamente il 'trattamento' e i 'blocchi'. Ciascuna cella sarà a sua volta riempita con il terzo fattore, denominato 'fattore principale'. Tutti e 3 i criteri devono avere lo stesso numero di livelli.

Per organizzare un disegno sperimentale di questo genere, ci risulta utile il software R già per organizzare il quadrato latino, in quanto non vi devono essere ripetizioni nell'ambito delle righe e delle colonne. Per ottenere un quadrato latino da utilizzare, usiamo la library(crossdes):


library(crossdes)
LatSqua <- des.MOLS(5,5)
LatSqua[1:5,1:5]

[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 2 3 4 5 1
[3,] 3 4 5 1 2
[4,] 4 5 1 2 3
[5,] 5 1 2 3 4


Abbiamo così ottenuto un possibile schema di quadrato latino da utilizzare. Una volta raccolti i dati, potremo passare alla fase analitica.

Supponiamo ad esempio di avere un appezzamento di terreno, e di dividerlo in 5 strisce e 5 colonne, così da ottenere 25 quadrati. Per ogni striscia viene utilizzata una profondità di aratura differente; in ogni colonna viene utilizzato un concime differente; per ogni combinazione viene utilizzato una diversa varietà di seme. La nostra domanda è: la produttività dipende dal seme, dalla profondità dell'aratura o dal tipo di concime usato? O da tutt'e 3 i fattori?

Adesso supponiamo di avere ottenuto il seguente quadrato latino:


trattamentoA trattamentoB trattamentoC trattamentoD trattamentoE
concime1 "A42" "C47" "B55" "D51" "E44"
concime2 "E45" "B54" "C52" "A44" "D50"
concime3 "C41" "A46" "D57" "E47" "B48"
concime4 "B56" "D52" "E49" "C50" "A43"
concime5 "D47" "E49" "A45" "B54" "C46"


I tre fattori sono quindi: tipo di concime (concime1:5), tipo di trattamento del terreno (trattamentoA:E), tipo di seme (A:E). I numeri rappresentano la produttività, in quintali, di un certo frutto.

Per analizzare questi dati in R dobbiamo creare un dataframe. Utilizziamo il seguente codice:


concime <- c(rep("concime1",1), rep("concime2",1), rep("concime3",1), rep("concime4",1), rep("concime5",1))
trattamento <- c(rep("tratA",5), rep("tratB",5), rep("tratC",5), rep("tratD",5), rep("tratE",5))
varieta <- c("A","E","C","B","D", "C","B","A","D","E", "B","C","D","E","A", "D","A","E","C","B", "E","D","B","A","C")
freq <- c(42,45,41,56,47, 47,54,46,52,49, 55,52,57,49,45, 51,44,47,50,54, 44,50,48,43,46)

mydata <- data.frame(trattamento, concime, varieta, freq)

mydata

trattamento concime varieta freq
1 tratA concime1 A 42
2 tratA concime2 E 45
3 tratA concime3 C 41
4 tratA concime4 B 56
5 tratA concime5 D 47
6 tratB concime1 C 47
7 tratB concime2 B 54
8 tratB concime3 A 46
9 tratB concime4 D 52
10 tratB concime5 E 49
11 tratC concime1 B 55
12 tratC concime2 C 52
13 tratC concime3 D 57
14 tratC concime4 E 49
15 tratC concime5 A 45
16 tratD concime1 D 51
17 tratD concime2 A 44
18 tratD concime3 E 47
19 tratD concime4 C 50
20 tratD concime5 B 54
21 tratE concime1 E 44
22 tratE concime2 D 50
23 tratE concime3 B 48
24 tratE concime4 A 43
25 tratE concime5 C 46


Per conferma, creiamo due matrici, con le varietà, e con le frequenze:


matrix(mydata$varieta, 5,5)

[,1] [,2] [,3] [,4] [,5]
[1,] "A" "C" "B" "D" "E"
[2,] "E" "B" "C" "A" "D"
[3,] "C" "A" "D" "E" "B"
[4,] "B" "D" "E" "C" "A"
[5,] "D" "E" "A" "B" "C"

matrix(mydata$freq, 5,5)

[,1] [,2] [,3] [,4] [,5]
[1,] 42 47 55 51 44
[2,] 45 54 52 44 50
[3,] 41 46 57 47 48
[4,] 56 52 49 50 43
[5,] 47 49 45 54 46


Osservate che queste due tabelle combinate, danno la tabella originale da cui abbiamo preso i dati, per cui abbiamo inserito correttamente i valori.

Prima di procedere con l'analisi della varianza di questo disegno a quadrato latino, è opportuno eseguire un boxplot, per avere un'idea del risultato che ci aspettiamo:


par(mfrow=c(2,2))
plot(freq ~ concime+trattamento+varieta, mydata)




Osservate che le medie sono realmente differenti considerando il fattore varieta, poco differenti considerando il fattore trattamento, praticamente uguali considerando il tipo di concime usato.
Verifichiamo queste osservazioni grafiche, con la tabella ANOVA:


myfit <- lm(freq ~ concime+trattamento+varieta, mydata)
anova(myfit)

Analysis of Variance Table

Response: freq
Df Sum Sq Mean Sq F value Pr(>F)
concime 4 17.760 4.440 0.7967 0.549839
trattamento 4 109.360 27.340 4.9055 0.014105 *
varieta 4 286.160 71.540 12.8361 0.000271 ***
Residuals 12 66.880 5.573
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Osservate la significatività (il p-value) dei valori F di Fisher per i 3 test che sono stati effettuati. Considerando il fattore concime, non c'è differenza tra le medie; considerando il fattore trattamento la significatività è moderata (p-value compreso tra 0.05 e 0.1), mentre considerando il fattore varietà, la differenza è altamente significativa.

Adesso procediamo con la diagnostica del modello, osservando la distribuzione dei residui:



La distribuzione dei residui del primo grafico è casuale; il Q-Q plot dei residui è distribuito sufficientemente sulla bisettrice, quindi concludo che i residui seguono una distribuzione normale, assunzione indispensabile per un modello lineare come il nostro.