Cu experimentele de teren, studiul corelației dintre trăsăturile observate poate să nu fie o sarcină ușoară. De exemplu, putem lua în considerare un experiment de genotip, desfășurat în blocuri complete randomizate, cu 27 de genotipuri de grâu și trei replici, în care au fost înregistrate mai multe trăsături, inclusiv randamentul (Randamentul) și greutatea unei mii de boabe (TKW). Am putea fi interesați să studiem corelația dintre aceste două trăsături, dar ar trebui să ne confruntăm cu două probleme fundamentale:
- conceptul de corelare într-o astfel de situație nu este unic, deoarece am putea lua în considerare fie corelația dintre măsurătorile plotului, fie corelația dintre reziduuri sau corelația dintre mediile genotipului sau corelația dintre mediile bloc;
- unitățile experimentale nu sunt independente, dar sunt grupate după genotip și bloc, ceea ce invalidează toate inferențe bazate pe ipoteza independenței.
Am tratat aceste două probleme (în special prima) într-un post anterior, unde am dat o soluție bazată pe metode tradiționale de analiză a datelor.
În această postare, aș dori să prezint o soluție mai avansată, care a fost susținută de Hans-Peter Piepho într-un manuscris relativ recent (Piepho, 2018). O astfel de soluție se bazează pe modele mixte și a fost implementată în SAS, prin utilizarea PROC MIXED. Am petrecut câteva ore „transportând” acele modele în R, ceea ce s-a dovedit a fi o sarcină grea, deși, până la urmă, se pare că am ajuns la o soluție acceptabilă, pe care aș vrea să o împărtășesc aici.
În primul rând, putem încărca setul de date „WheatQuality”, care este disponibil în pachetul „statforbiology”; constă din 81 de înregistrări (parcele) a 6 variabile, adică factorii Genotip și Block, precum și cele patru răspunsuri înălțime, TKW, greutate pe hectolitru și randament. Codul de mai jos încarcă pachetele necesare, datele și transformă variabila numerică „Block” într-un factor.
rm(list = ls())
library(statforbiology)
library(dplyr)
library(tidyr)
library(sommer)
library(nlme)
#
# Loading data
dataset <- getAgroData("WheatQuality") |>
mutate(Block = factor(Block),
Genotype = factor(Genotype))
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
Piepho (2018) a arătat că, pentru un experiment ca acesta, toți coeficienții de corelație pot fi estimați prin codificarea unui model mixt cu răspunsuri multiple, după cum urmează:
( Y_{ijk} = mu_i + beta_{ik} + tau_{ij} + epsilon_{ijk})
unde (Y_{ijk}) este răspunsul pentru trăsătură (i)genotipul (j) si blocul (k), (mu_i) este media pentru trăsătură (i), (beta_{ik}) este efectul blocului (k) și trăsătură (i), (tau_{ij}) este efectul genotipului (j) pentru trăsătură (i) şi (epsilon_{ijk}) este rezidual pentru trăsătură (i)genotipul (j) si blocul (k).
În modelul de mai sus, reziduurile (epsilon_{ijk}) trebuie să fie distribuite în mod normal și heteroscedastice, cu variații specifice trăsăturilor. În plus, reziduurile aparținând aceluiași diagramă (cele două valori observate pentru cele două trăsături) trebuie corelate (corelarea erorilor).
Hans-Peter Piepho, în lucrarea sa, a prezentat ideea că efectele „genotip” și „bloc” pentru cele două trăsături pot fi luate ca aleatoare, ceea ce are sens, deoarece, pentru această aplicație, ne interesează în principal varianțele și covarianțele. Ambele efecte aleatoare (pentru genotip și pentru termenii bloc) trebuie să fie heteroscedastice (componente de varianță specifice trăsăturii) și trebuie să existe o corelație între cele două trăsături.
Trebuie remarcat că, pentru alte aplicații, genotipul și efectele de bloc (în special cele dintâi) ar putea fi considerate mai bine ca fiind fixe, dar nu vom urmări o astfel de idee în acest post.
Din câte știu, nu există nicio modalitate de a potrivi un model atât de complex cu pachetul „nlme” și cu funcția „lme()” aferentă (voi da un indiciu mai târziu, pentru un model mai simplu). Într-o postare anterioară la acest link, am oferit o soluție bazată pe pachetul „asreml” (Butler et al., 2018), dar aceasta este o opțiune plătită. În vremuri mai recente am descoperit pachetul „sommer” (Covarrubias-Pazaran, 2016), care pare a fi foarte util și este potrivit pentru a face față cerințelor acestei postări. Funcția cheie a lui „sommer” este mmer()și, pentru a se potrivi cu modelul de mai sus, trebuie să specificăm următoarele componente.
- Variabilele de răspuns. Când setăm un model bivariat cu „sommer”, putem „cbind()” Randament și TKW.
- Modelul fix, care nu conține niciun efect, ci interceptarea (în mod implicit, mijloacele pentru cele două efecte sunt estimate separat, ca în Piepho, 2018).
- Modelul aleatoriu, care este compus din efectele „genotip” și „bloc”. Pentru ambele, am specificat o matrice de covarianță generală nestructurată, astfel încât să putem estima două componente de varianță diferite (una per trăsătură) și o componentă de covarianță. Codarea rezultată este „~ vsr(usr(Genotip)) + vsr(usr(Block))”.
- Structura reziduală, în care cele două observații din aceeași diagramă sunt heteroscedastice și corelate. Această structură este instalată implicit și nu necesită nicio codificare specifică.
Apelul model este:
mod.bimix <- mmer(cbind(Yield, TKW) ~ 1,
random = ~ vsr(usr(Genotype)) + vsr(usr(Block)),
rcov = ~ vsr(units),
data = dataset,
verbose = FALSE, dateWarning = FALSE)
bimix.obj <- summary(mod.bimix)
coefs <- bimix.obj$varcomp
coefs
## VarComp VarCompSE Zratio Constraint
## u:Genotype.Yield-Yield 77.6342608 22.0978545 3.5132035 Positive
## u:Genotype.Yield-TKW 38.8320973 15.0930429 2.5728475 Unconstr
## u:Genotype.TKW-TKW 53.8613303 15.3539585 3.5079768 Positive
## u:Block.Yield-Yield 3.7104682 3.9363372 0.9426195 Positive
## u:Block.Yield-TKW -0.2428322 1.9074202 -0.1273092 Unconstr
## u:Block.TKW-TKW 1.6684549 1.8343512 0.9095613 Positive
## u:units.Yield-Yield 6.0939217 1.1951482 5.0988836 Positive
## u:units.Yield-TKW 0.1635821 0.7242898 0.2258518 Unconstr
## u:units.TKW-TKW 4.4718011 0.8770118 5.0989065 Positive
Caseta de mai sus arată rezultatele despre parametrii varianță-covarianță. Pentru a obține corelațiile, am folosit metoda delta, așa cum este implementată în gnlht() funcția pachetului „statforbiology” (pachetul însoțitor pentru acest blog). În primul rând am extras parametrii de varianță împreună cu matricea de covarianță pentru parametrii de varianță din obiectul model mixt. Pentru simplitate, am atribuit nume simple coeficienților (V1, V2, … Vn), în funcție de ordonarea lor în ieșirea modelului.
# Correlation between genotype means
coefsVec <- coefs(,1)
vcovMat <- mod.bimix$sigmaSE # Variance-covariance for varcomp
names(coefsVec) <- paste("V", 1:9, sep = "")
gnlht(coefsVec, func = list(~ V2 / (sqrt(V1)*sqrt(V3) )),
vcov. = as.matrix(vcovMat),
parameterNames = paste("V", 1:9, sep = ""))
## Form Estimate SE Z-value p-value
## 1 V2/(sqrt(V1) * sqrt(V3)) 0.6005174 0.1306699 4.595684 4.313326e-06
#
# Correlation between block means
gnlht(coefsVec, func = list(~ V5 / (sqrt(V4)*sqrt(V6) ) ),
vcov. = as.matrix(vcovMat),
parameterNames = paste("V", 1:9, sep = ""))
## Form Estimate SE Z-value p-value
## 1 V5/(sqrt(V4) * sqrt(V6)) -0.09759658 0.7571256 0.1289041 0.8974335
#
# Correlation of residuals
gnlht(coefsVec, func = list(~ V8 / (sqrt(V7)*sqrt(V9) )),
vcov. = as.matrix(vcovMat),
parameterNames = paste("V", 1:9, sep = ""))
## Form Estimate SE Z-value p-value
## 1 V8/(sqrt(V7) * sqrt(V9)) 0.03133619 0.1385421 0.2261854 0.8210572
Vedem că estimările sunt foarte apropiate de cele obținute prin utilizarea coeficienților de corelație ai lui Pearson (vezi postarea mea anterioară). Avantajul acestei soluții de model mixt este că putem testa și ipoteze într-un mod relativ fiabil. De exemplu, putem analiza testele Wald din rezultatul de mai sus pentru a judeca semnificația corelațiilor și a concluziona că numai mediile genotipului sunt corelate semnificativ între ele.
Având în vedere că mediile bloc nu sunt corelate, dacă am fi dispuși să luăm efectul „bloc” ca fix, am putea potrivi un model mixt bivariat și cu pachetul „nlme” și cu funcția lme() (Pinheiro și Bates, 2000). Cu toate acestea, ar trebui să aruncăm modelul ca un model „univariat”.
În acest scop, cele două variabile (Randament și TKW) trebuie să fie îngrămădite și trebuie introdus un nou factor („trăsătură”) pentru a identifica observațiile pentru cele două trăsături. Un alt factor este de asemenea necesar pentru a identifica diferitele parcele, adică unitățile de observație. Pentru a efectua o astfel de restructurare, am folosit pivot_longer() funcția în pachetul „tidyr” (Wickham et al., 2024) și a atribuit numele „Y” variabilei de răspuns, care este acum compusă din cele două variabile originale Yield și TKW, una peste alta.
dataset$Plot <- 1:81
mdataset <- dataset |>
select(-Whectol, -Height) |>
pivot_longer(names_to = "Trait", values_to = "Y", cols = c("Yield", "TKW")) |>
mutate(Trait = factor(Trait))
head(mdataset)
## # A tibble: 6 × 5
## Genotype Block Plot Trait Y
##
## 1 arcobaleno 1 1 Yield 64.4
## 2 arcobaleno 1 1 TKW 44.5
## 3 arcobaleno 2 2 Yield 60.6
## 4 arcobaleno 2 2 TKW 42.8
## 5 arcobaleno 3 3 Yield 59.4
## 6 arcobaleno 3 3 TKW 42.7
tail(mdataset)
## # A tibble: 6 × 5
## Genotype Block Plot Trait Y
##
## 1 vitromax 1 79 Yield 54.4
## 2 vitromax 1 79 TKW 41.6
## 3 vitromax 2 80 Yield 51.0
## 4 vitromax 2 80 TKW 43.6
## 5 vitromax 3 81 Yield 48.8
## 6 vitromax 3 81 TKW 43.1
Modelul fix este:
Y ~ Trait*Block
Pentru a introduce o matrice varianță-covarianță total nestructurată pentru efectul aleator, am folosit constructul „pdMat”:
random = list(Genotype = pdSymm(~Trait - 1))
În ceea ce privește reziduurile, heteroscedasticitatea poate fi inclusă folosind argumentul „weight()” și funcția de varianță „varIdent”, care permite o abatere standard diferită per trăsătură:
weight = varIdent(form = ~1|Trait)
De asemenea, am permis corelarea reziduurilor în cadrul fiecărei parcele, utilizând argumentul „corelație” și specificând o formă generală de corelare „corSymm()”. Parcele sunt imbricate în genotipuri, așa că am folosit un operator de cuibărit, după cum urmează:
correlation = corSymm(form = ~1|Genotype/Plot)
Apelul model final este:
mod <- lme(Y ~ Trait*Block,
random = list(Genotype = pdSymm(~Trait - 1)),
weight = varIdent(form=~1|Trait),
correlation = corCompSymm(form=~1|Genotype/Plot),
data = mdataset)
Recuperarea matricelor varianță-covarianță necesită un efort, deoarece funcția „getVarCov()” nu pare să funcționeze corect cu acest model cu mai multe straturi. În primul rând, putem prelua matricea varianță-covarianță pentru efectul aleator al genotipului (G) și matricea de corelație corespunzătoare.
#Variance structure for random effects obj <- mod G <- matrix( as.numeric(getVarCov(obj)), 2, 2 ) G ## (,1) (,2) ## (1,) 53.86053 38.83124 ## (2,) 38.83124 77.63402 cov2cor(G) ## (,1) (,2) ## (1,) 1.0000000 0.6005096 ## (2,) 0.6005096 1.0000000
În al doilea rând, putem prelua matricea „condițională” varianță-covarianță (R), care descrie corelarea erorilor:
#Conditional variance-covariance matrix (residual error) V <- corMatrix(obj$modelStruct$corStruct)((1)) #Correlation for residuals sds <- 1/varWeights(obj$modelStruct$varStruct)(1:2) sds <- obj$sigma * sds #Standard deviations for residuals (one per trait) R <- t(V * sds) * sds #Going from correlation to covariance R ## (,1) (,2) ## (1,) 6.0939152 0.1634968 ## (2,) 0.1634968 4.4718077 cov2cor(R) ## (,1) (,2) ## (1,) 1.00000000 0.03131984 ## (2,) 0.03131984 1.00000000
Matricea de corelație totală se obține pur și simplu ca sumă a lui G și R:
Tr <- G + R cov2cor(Tr) ## (,1) (,2) ## (1,) 1.000000 0.555787 ## (2,) 0.555787 1.000000
Vedem că aceleași rezultate pot fi obținute folosind „sommer” și considerând efectul de blocare ca fiind fix, deși codarea de mai jos este mult mai îngrijită!
mod.bimix5 <- mmer(cbind(Yield, TKW) ~ Block,
random = ~ vsr(usr(Genotype)),
data = dataset,
verbose = FALSE, dateWarning = FALSE)
mod.bimix5$sigma
## $`u:Genotype`
## Yield TKW
## Yield 77.63425 38.83209
## TKW 38.83209 53.86133
##
## $units
## Yield TKW
## Yield 6.0939217 0.1635824
## TKW 0.1635824 4.4718011
Sper că a fost util… dacă aveți soluții mai bune, aș fi bucuros să le învăț; Vă rog, trimiteți-mi un rând la adresa de mai jos. Mulțumim pentru lectură și codificare fericită!
Și… nu uitați să verificați noua mea carte!
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
