(Acest articol a fost publicat pentru prima dată pe Statforbiologieși cu amabilitate a contribuit la R-bloggeri). (Puteți raporta probleme legate de conținutul acestei pagini aici)
Doriți să vă distribuiți conținutul pe R-bloggeri? dați clic aici dacă aveți un blog, sau aici dacă nu aveți.
Unul dintre scopurile fundamentale ale experimentelor de teren este explorarea relațiilor dintre mai multe trăsături experimentale, care au fost înregistrate/măsurate în totalitate pentru toate unitățile experimentale. De exemplu, gândiți-vă la un experiment de genotip cu 27 de genotipuri de grâu de iarnă într-un proiect de bloc complet randomizat cu trei replici. Dacă ne propunem să apreciem randamentul și calitatea randamentului, trebuie să măsurăm, în fiecare parcelă, mai multe trăsături, cum ar fi randamentul cerealelor și greutatea unei mii de boabe (TKW). Pe lângă analiza separată a acestor două variabile, am putea fi interesați să descriem relațiile lor reciproce, folosind coeficientul de corelație al lui Pearson. Să facem asta, utilizând setul de date „WheatQuality”, care este disponibil în pachetul „statforbiology”; avem 81 de parcele și la fel de multe cupluri de măsuri în total.
library(Hmisc)
library(dplyr)
library(statforbiology)
dataset <- getAgroData("WheatQuality")
head(dataset)
## Genotype Block Height TKW Whectol Yield
## 1 arcobaleno 1 90 44.5 83.2 64.40
## 2 arcobaleno 2 90 42.8 82.2 60.58
## 3 arcobaleno 3 88 42.7 83.1 59.42
## 4 baio 1 80 40.6 81.8 51.93
## 5 baio 2 75 42.7 81.3 51.34
## 6 baio 3 76 41.1 81.1 47.78
Acum, luăm în considerare întregul lot de măsuri și calculăm coeficientul de corelație dintre randament și TKW; folosim și un test t pentru ipoteza că coeficientul de corelație este egal cu 0. Acesta este rezultatul:
ra <- with(dataset, rcorr(Yield, TKW) ) # # Correlation coefficient ra$r(1,2) ## (1) 0.540957 # # P-level ra$P(1,2) ## (1) 1.850931e-07
Concluzia este: coeficientul de corelație este pozitiv și semnificativ diferit de 0, ceea ce susține ideea că parcelele cu o valoare mare a randamentului tind să aibă și o valoare mare a TKW.
O astfel de abordare se vede foarte des și în reviste de rang înalt, dar: este corect? Putem avea încredere în concluzii? Răspunsul este nu, nu complet. Și iată de ce. Vă rugăm să rețineți că în această postare urmăresc punctul susținut de Bland și Altman (1994), folosind un exemplu care este mai strâns legat de știința plantelor/culturii.
Coeficientul de corelație este o măsură a variabilității comune și face câteva ipoteze de bază (de exemplu, liniaritatea), pe care le-am detaliat în altă parte (vezi aici). În special, testul t necesită ca toate cuplurile de observații să fie luate pe subiecți independenți.
În cazul nostru, cele 81 de parcele nu sunt independente, deoarece sunt grupate pe bloc și genotip, astfel încât observațiile din același bloc sau pentru același genotip sunt mai asemănătoare decât observațiile din blocuri diferite și genotipuri diferite (corelație intraclasă). În consecință, testul t este invalid și nu avem voie să concluzionăm că corelația este semnificativă, deoarece o astfel de concluzie nu este neapărat susținută de date.
Din păcate, o astfel de lipsă de independență se referă la toate experimentele proiectate, ceea ce pare a fi un fapt care trece în mare parte neobservat. În general, coeficienții de corelație calculați pentru datele brute din experimentele proiectate nu pot fi supuși testelor t, fără a lua în considerare efectele tuturor factorilor de tratament și grupare.
În afară de testarea statistică, coeficientul de corelație pe care l-am calculat mai sus (0,541) își menține valoarea ca măsură descriptivă a variabilității articulare. Dar suntem siguri că utilizarea datelor brute observate este cea mai bună cale de a merge?
În primul rând, un astfel de coeficient de corelație evidențiază că randamentul și TKW sunt corelate pozitiv, dar este acest lucru adevărat pentru toate genotipurile și pentru toate blocurile? Să verificăm asta.
# Correlations within genotypes wgCors <- by(dataset, list(dataset$Genotype), function(x) cor(x$Yield, x$TKW)) wgCors ## : arcobaleno ## (1) 0.9847228 ## ------------------------------------------------------------ ## : baio ## (1) 0.1611952 ## ------------------------------------------------------------ ## : claudio ## (1) -0.9993872 ## ------------------------------------------------------------ ## : colorado ## (1) 0.9837293 ## ------------------------------------------------------------ ## : colosseo ## (1) 0.4564855 ## ------------------------------------------------------------ ## : creso ## (1) -0.5910193 ## ------------------------------------------------------------ ## : duilio ## (1) -0.988233 ## ------------------------------------------------------------ ## : gianni ## (1) -0.7603802 ## ------------------------------------------------------------ ## : giotto ## (1) 0.9520211 ## ------------------------------------------------------------ ## : grazia ## (1) 0.4980828 ## ------------------------------------------------------------ ## : iride ## (1) 0.7563338 ## ------------------------------------------------------------ ## : meridano ## (1) 0.1174342 ## ------------------------------------------------------------ ## : neodur ## (1) 0.4805871 ## ------------------------------------------------------------ ## : orobel ## (1) 0.8907754 ## ------------------------------------------------------------ ## : pietrafitta ## (1) -0.9633891 ## ------------------------------------------------------------ ## : portobello ## (1) 0.9210135 ## ------------------------------------------------------------ ## : portofino ## (1) -0.9900764 ## ------------------------------------------------------------ ## : portorico ## (1) 0.1394211 ## ------------------------------------------------------------ ## : preco ## (1) 0.9007067 ## ------------------------------------------------------------ ## : quadrato ## (1) -0.5840238 ## ------------------------------------------------------------ ## : sancarlo ## (1) -0.646067 ## ------------------------------------------------------------ ## : simeto ## (1) -0.4051779 ## ------------------------------------------------------------ ## : solex ## (1) -0.6066363 ## ------------------------------------------------------------ ## : terrabianca ## (1) -0.4076416 ## ------------------------------------------------------------ ## : verdi ## (1) 0.5801404 ## ------------------------------------------------------------ ## : vesuvio ## (1) -0.7797493 ## ------------------------------------------------------------ ## : vitromax ## (1) -0.8056514 # Correlations within blocks wbCors <- by(dataset, list(dataset$Block), function(x) cor(x$Yield, x$TKW)) wbCors ## : 1 ## (1) 0.5998137 ## ------------------------------------------------------------ ## : 2 ## (1) 0.539999 ## ------------------------------------------------------------ ## : 3 ## (1) 0.5370398
În ceea ce privește genotipurile, avem 27 de coeficienți de corelație, variind de la -0,999 la 0,985. Avem doar trei cupluri de măsurători per genotip și este clar că acest lucru nu este mult, pentru a estima fiabil un coeficient de corelație. Cu toate acestea, este suficient să fii suspicios cu privire la amploarea corelației dintre randament și TKW, deoarece poate depinde de genotip.
Pe de altă parte, corelația în cadrul blocurilor este mai constantă, independentă de bloc și similară cu corelația dintre parcele.
Poate fi interesant să obținem o estimare a corelației medii în cadrul grupului. În acest scop, putem efectua două ANOVA separate (unul pe trăsătură), incluzând toate efectele relevante (blocuri și genotipuri) și să calculăm corelația dintre reziduuri:
mod1 <- lm(Yield ~ factor(Block) + Genotype, data = dataset) mod2 <- lm(TKW ~ factor(Block) + Genotype, data = dataset) wCor <- rcorr(residuals(mod1), residuals(mod2)) wCor$r ## x y ## x 1.00000000 0.03133109 ## y 0.03133109 1.00000000 wCor$P ## x y ## x NA 0.7812693 ## y 0.7812693 NA
Corelația medie în cadrul grupului este foarte mică și nesemnificativă. Să ne gândim la această corelație: reziduurile reprezintă valorile de randament și TKW pentru toate parcelele, odată ce efectele blocurilor și genotipurilor au fost eliminate. O corelație ridicată a reziduurilor ar însemna că, lăsând deoparte efectele blocului și genotipului căruia îi aparține, o parcelă cu o valoare mare asupra randamentului prezintă și o valoare mare a TKW. Concluzia corectă, într-adevăr, este că, odată ce efectele tratamentului și factorii de grupare au fost eliminate, nu există nicio corelație între trăsăturile studiate.
Dar, acest lucru nu este suficient pentru a opri investigațiile noastre. De asemenea, am putea lua în considerare mijloacele pentru genotipuri/blocuri și am putea vedea dacă acestea sunt corelate. Blocurile și genotipurile sunt independente și, în principiu, testarea semnificației este permisă. Cu toate acestea, acest lucru nu este recomandat cu mijloace de blocare, deoarece trei date sunt prea puține pentru a face teste.
means <- aggregate(dataset(,c(4,6)), list(dataset$Genotype), mean) rgPear <- rcorr( as.matrix(means(,2:3)) ) rgPear$r ## TKW Yield ## TKW 1.0000000 0.5855966 ## Yield 0.5855966 1.0000000 rgPear$P ## TKW Yield ## TKW NA 0.00133149 ## Yield 0.00133149 NA means <- aggregate(dataset(,c(4,6)), list(dataset$Block), mean) rbPear <- cor( as.matrix(means(,2:3)) ) rbPear ## TKW Yield ## TKW 1.00000000 -0.08812544 ## Yield -0.08812544 1.00000000
Observăm că corelația dintre mediile genotipului este ridicată și semnificativă. Dimpotrivă, corelația dintre mediile blocului este aproape de 0.
În această etapă, s-ar putea să fii confuz… Să încercăm să curățăm ceața.
O analiză naivă a corelației asupra datelor brute ne-a adus la concluzia că randamentul și TKW sunt corelate semnificativ. Din păcate, o astfel de concluzie este prea simplistă și nu este susținută cu greu de date.
Un proces mai atent de analiză a datelor nu trebuie să neglijeze faptul că în fiecare experiment de teren proiectat avem grupuri de subiecți și există mai multe tipuri posibile de corelație:
- corelația dintre măsurătorile grafice (corelație naivă)
- corelația dintre reziduuri
- corelația dintre mijloacele de tratament/blocare
Toate aceste corelații ar trebui investigate și utilizate pentru interpretare.
- Corelația „naivă” dintre măsurătorile parcelei este foarte ușor de calculat, dar este extrem de înșelătoare. Într-adevăr, nu ține cont de efectele de tratament/blocare și nu permite testarea ipotezelor, subiecții nefiind independenți. Cred că este rezonabil să sugerăm că coeficientul de corelație „naiv” nu este niciodată utilizat pentru interpretare.
- Corelația dintre reziduuri este o măsură sigură a variației articulațiilor, deoarece designul experimental este referit corespunzător, prin eliminarea efectelor tratamentelor/blocurilor. În acest exemplu, reziduurile nu sunt corelate. Analize ulterioare arată că corelația dintre randament și TKW, dacă există, poate depinde de genotip, în timp ce nu depinde de bloc.
- Corelația dintre mijloacele tratament/blocare permite să se evalueze dacă efectele tratament/blocare asupra celor două trăsături sunt corelate. În acest caz, deși nu avem voie să concluzionăm că randamentul și TKW sunt, în general, corelate, putem concluziona că genotipurile cu un nivel ridicat de randament prezintă și un nivel ridicat de TKW.
Piepho (2018), în cadrul modelelor mixte, a propus o metodă elegantă și avansată de abordare a analizelor de corelare a datelor din experimentele de teren proiectate. Dar asta merită încă o postare.
Mulțumesc pentru lectură și codificare fericită!
(PS Dacă v-a plăcut această postare, consultați cartea mea „Metode de cercetare pe teren în agricultură: o introducere cu R” prin linkul de mai jos!)
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
