Cu experimentele de teren, studierea 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, stabilit în blocuri complete randomizate, cu 27 de genotipuri de grâu și trei replici, unde au fost înregistrate mai multe trăsături, inclusiv randamentul (randamentul) și greutatea a mii de nuclee (TKW). Am putea fi interesați să studiem corelația dintre cele două trăsături, dar ar trebui să ne confruntăm cu două probleme fundamentale:
- Conceptul de corelație într -o astfel de setare nu este unic, deoarece am putea lua în considerare corelația dintre măsurătorile parcelei, fie corelația dintre reziduuri sau corelația dintre mijloacele genotipului sau corelația dintre mijloacele blocului;
- Unitățile experimentale nu sunt independente, dar sunt grupate de genotip și bloc, care invalidează toate inferențele pe baza presupunerii independenței.
M -am ocupat de aceste două probleme (în special de prima) într -o postare anterioară, unde am dat o soluție bazată pe metodele tradiționale de analize de date.
Î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, care s -au dovedit a fi o sarcină dificilă, deși, până la urmă, par să fi ajuns la o soluție acceptabilă, pe care aș dori să o împărtășesc aici.
În primul rând, putem încărca setul de date „WheatQuality.csv”, care este disponibil este disponibil pe un depozit de internet; Este format din 81 de înregistrări (parcele) de 6 variabile, adică factorii de genotip și bloc, 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ă „bloc” într -un factor.
rm(list = ls())
library(dplyr)
library(sommer)
library(nlme)
#
# Loading data
dataset <- read.csv("https://www.casaonofri.it/_datasets/WheatQuality.csv") |>
mutate(Block = factor(Block))
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 estimate prin codificarea unui model mixt cu mai multe răspunsuri, 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 ) și 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 reziduul pentru trăsătură (i )genotipul (j ) și blocul (k ).
În modelul de mai sus, reziduurile ( epsilon_ {ijk} ) trebuie să fie în mod normal distribuit și heteroscedastic, cu variații specifice trăsăturilor. Mai mult, reziduurile aparținând aceluiași complot (cele două valori observate pentru cele două trăsături) trebuie să fie corelate (corelarea erorilor).
Hans-Peter Piepho, în lucrarea sa, a prezentat ideea că efectele „genotipului” și „bloc” pentru cele două trăsături pot fi luate ca fiind aleatorii, ceea ce are sens, deoarece, pentru această aplicație, suntem interesați în principal de variații și covarianțe. Ambele efecte aleatorii (pentru genotip și pentru termenii bloc) trebuie să fie heteroscedastice (componente de varianță specifică trăsăturilor) și trebuie să existe o corelație între cele două trăsături.
Trebuie menționat că, pentru alte aplicații, efectele genotipului și blocului (în special primele) ar putea fi mai bine considerate fixate, dar nu vom urmări o astfel de idee în acest post.
În conformitate cu cunoștințele mele, nu există nicio modalitate de a se potrivi unui model atât de complex cu pachetul „NLME” și funcția „LME ()” înrudită (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 și colab., 2018), dar aceasta este o opțiune plătită. În vremurile mai recente am descoperit pachetul „Sommer” (Covarrubias-Pazaran, 2016), care pare a fi foarte util și este potrivit pentru a face față cerințelor acestui post. Funcția cheie a „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 stabilim un model bivariat cu „Sommer”, putem „cbind ()” randament și TKW.
- Modelul fix, care nu conține niciun efecte, dar 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 generală de covarianță de varianță nestructurată, astfel încât să putem estima două componente de varianță diferite (una pe trăsătură) și o componentă de covarianță. Codificarea rezultată este „~ vsr (usr (genotip)) + vsr (usr (bloc))”.
- Structura reziduală, unde cele două observații din același complot sunt heteroscedestice și corelate. Această structură este montată în mod implicit și nu necesită nicio codificare specifică.
Apelul model este:
mod.bimix <- mmer(cbind(Yield, TKW) ~ 1,
random = ~ vsr(usr(Genotype)) + vsr(usr(Block)),
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
## units.Yield-Yield 6.0939217 1.1951482 5.0988836 Positive
## units.Yield-TKW 0.1635821 0.7242898 0.2258518 Unconstr
## units.TKW-TKW 4.4718011 0.8770118 5.0989065 Positive
Caseta de mai sus arată rezultatele despre parametrii variației-covarianță. Pentru a obține corelațiile, am folosit funcția vpredict() și a specificat combinațiile necesare parametrilor de varianță-covarianță. Este necesar să ne amintim că estimările, în „Sommer”, sunt numite V1, V2,… VN, în funcție de ordonarea lor în ieșirea modelului.
# Correlation between genotype means vpredict(mod.bimix, rg ~ V2 / (sqrt(V1)*sqrt(V3) )) ## Estimate SE ## rg 0.6005174 0.1306699 # Correlation between block means vpredict(mod.bimix, rt ~ V5 / (sqrt(V4)*sqrt(V6) )) ## Estimate SE ## rt -0.09759658 0.7571256 vpredict(mod.bimix, rt ~ V8 / (sqrt(V7)*sqrt(V9) )) ## Estimate SE ## rt 0.03133619 0.1385421
Vedem că estimările sunt foarte apropiate de cele obținute prin utilizarea coeficienților de corelație 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, am testat ipoteza că reziduurile nu sunt corelate de:
- codificarea unui model redus în care reziduurile sunt heteroscedastice și independente și
- Compararea acestui model redus cu modelul complet, printr-un test de raport de probabilitate bazat pe Reml (LRT).
Înlăturarea corelației reziduurilor se face cu ușurință, folosind codul de mai jos.
mod.bimix2 <- mmer(cbind(Yield, TKW) ~ 1,
random = ~ vsr(usr(Genotype)) + vsr(usr(Block)),
rcov = ~ vsr(units, Gtc = diag(2)),
data = dataset,
verbose = FALSE, dateWarning = FALSE)
# Likelihood ratio test
LRT <- anova(mod.bimix, mod.bimix2)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 15 -62.41068 -56.23549 33.20534
## mod2 14 -62.35960 -56.18441 33.17980 0.05108 1 0.8212
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De asemenea, am încercat să reduc în continuare modelul pentru a testa semnificația corelației dintre mijloacele blocului și mijloacele genotipului.
# Removing correlation between block means
mod.bimix3 <- mmer(cbind(Yield, TKW) ~ 1,
random = ~ vsr(usr(Genotype)) + vsr(Block, Gtc = diag(2)),
rcov = ~ vsr(units, Gtc = diag(2)),
data = dataset,
verbose = FALSE, dateWarning = FALSE)
# Likelihood ratio test
LRT <- anova(mod.bimix, mod.bimix3)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 15 -62.41068 -56.23549 33.20534
## mod2 13 -62.34401 -56.16882 33.17200 0.06667 2 0.96722
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Removing correlation between genotype means
mod.bimix4 <- mmer(cbind(Yield, TKW) ~ 1,
random = ~ vsr(Genotype, Gtc = diag(2)) + vsr(Block, Gtc = diag(2)),
rcov = ~ vsr(units, Gtc = diag(2)),
data = dataset,
verbose = FALSE, dateWarning = FALSE)
# Likelihood ratio test
LRT <- anova(mod.bimix, mod.bimix4)
## Likelihood ratio test for mixed models
## ==============================================================
## Df AIC BIC loLik Chisq ChiDf PrChisq
## mod1 15 -62.41068 -56.23549 33.20534
## mod2 12 -51.42519 -45.25000 27.71260 10.98549 3 0.0118 *
## ==============================================================
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vedem că doar mijloacele genotipului sunt corelate semnificativ între ele.
Având în vedere că mijloacele blocului nu sunt corelate, dacă am fi dispuși să luăm efectul „bloc” așa cum s -a fixat, am putea să încadrăm un model mixt bivariat, de asemenea, cu pachetul „NLME” și funcția lme() (Pinheiro și Bates, 2000). Cu toate acestea, ar trebui să aruncăm modelul ca model univariat.
În acest scop, trebuie să fie acumulate cele două variabile (randament și TKW) și trebuie introdus un nou factor („trăsături”) 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 melt() Funcția în pachetul „ResHape2” (Wickham, 2007) și a atribuit numele „Y” variabilei de răspuns, care este acum compusă de cele două variabile originale Randament și TKW, unul pe celălalt.
dataset$Plot <- 1:81
mdataset <- reshape2::melt(dataset(,c(-3,-5)), variable.name= "Trait", value.name="Y", id=c("Genotype", "Block", "Plot"))
mdataset$Block <- factor(mdataset$Block)
head(mdataset)
## Genotype Block Plot Trait Y
## 1 arcobaleno 1 1 TKW 44.5
## 2 arcobaleno 2 2 TKW 42.8
## 3 arcobaleno 3 3 TKW 42.7
## 4 baio 1 4 TKW 40.6
## 5 baio 2 5 TKW 42.7
## 6 baio 3 6 TKW 41.1
tail(mdataset)
## Genotype Block Plot Trait Y
## 157 vesuvio 1 76 Yield 54.37
## 158 vesuvio 2 77 Yield 55.02
## 159 vesuvio 3 78 Yield 53.41
## 160 vitromax 1 79 Yield 54.39
## 161 vitromax 2 80 Yield 50.97
## 162 vitromax 3 81 Yield 48.83
Modelul fix este:
Y ~ Trait*Block
Pentru a introduce o matrice de varianță-covarianță de varianță total nestructurată pentru efectul aleatoriu, am folosit construcția „PDMAT”:
random = list(Genotype = pdSymm(~Trait - 1))
În legătură cu reziduurile, heteroscedasticitatea poate fi inclusă folosind argumentul „greutate ()” și funcția de varianță „variată”, care permite o abatere standard diferită pe trăsătură:
weight = varIdent(form = ~1|Trait)
De asemenea, am permis corelarea reziduurilor în cadrul fiecărei comploturi, folosind argumentul „corelație” și specificând o formă de corelație generală „Corsymm ()”. Parcelele sunt cuibărite în genotipuri, astfel am folosit un operator de cuibărit, după cum urmează:
correlation = corSymm(form = ~1|Genotype/Plot)
Apelul final al modelului 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 matricilor de varianță-covarianță are nevoie de un anumit efort, deoarece funcția „getVarcov ()” nu pare să funcționeze corect cu acest model multi-strat. În primul rând, putem prelua matricea de varianță-covarianță pentru efectul aleatoriu 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.86081 38.83246 ## (2,) 38.83246 77.63485 cov2cor(G) ## (,1) (,2) ## (1,) 1.0000000 0.6005237 ## (2,) 0.6005237 1.0000000
În al doilea rând, putem prelua matricea de covarianță de varianță „condițională” (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,) 4.4718234 0.1635375 ## (2,) 0.1635375 6.0939251 cov2cor(R) ## (,1) (,2) ## (1,) 1.00000000 0.03132756 ## (2,) 0.03132756 1.00000000
Matricea de corelație totală este obținută pur și simplu ca suma lui G și R:
Tr <- G + R cov2cor(Tr) ## (,1) (,2) ## (1,) 1.0000000 0.5579906 ## (2,) 0.5579906 1.0000000
Vedem că aceleași rezultate pot fi obținute folosind „Sommer” și în ceea ce privește efectul bloc, deși codificarea este mai jos este mult mai naută!
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ă acest lucru a fost util … dacă aveți soluții mai bune, aș fi fericit să le învăț; Vă rog, aruncați -mi o linie la adresa de mai jos. Vă mulțumim pentru lectură și codificare fericită!
Andrea Onofri
Departamentul de Științe Agricole, Alimentare și Mediu
Universitatea din Perugia (Italia)
(e -mail protejat)
