--- title: "Fitting genotype by environment models in sommer" author: "Giovanny Covarrubias-Pazaran" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{GxE models in sommer} %\VignetteEngine{knitr::knitr} %\VignetteEncoding{UTF-8} --- The sommer package was developed to provide R users with a flexible univariate and multivariate linear mixed-model solver. The current package provides two complementary REML formulations. The `mmes()` function uses Henderson's mixed model equations and a sparse Average Information REML algorithm: the mixed-model coefficient matrix is factorized with sparse LDLT methods, selected inverse elements required by REML calculations can be obtained with Takahashi sparse-inverse recursions, and the factorization is reused for mixed-model solutions and variance-parameter derivatives. The `mmer()` function uses the marginal-covariance MNR formulation, in which REML calculations are organized around the observation covariance matrix \(V\) and the REML projection matrix \(P\). The relative efficiency of the two formulations depends on model dimensions, sparsity, and covariance structure rather than only on whether \(p>n\) or \(n>p\). The core numerical algorithms are coded in C++ using Armadillo and Eigen. The package allows users to specify flexible variance-covariance structures for random and residual effects and to obtain quantities such as REML variance-covariance estimates, BLUPs, BLUEs, residuals, fitted values, and prediction error variance information. The purpose of this vignette is to show how to fit different genotype by environment (GxE) models using the sommer package: 1) Single environment model 2) Multienvironment model: Main effect model 3) Multienvironment model: Diagonal model (DG) 4) Multienvironment model: Compund symmetry model (CS) 5) Multienvironment model: Unstructured model (US) 6) Multienvironment model: Random regression model (RR) 7) Multienvironment model: Other covariance structures for GxE 8) Multienvironment model: Finlay-Wilkinson regression 9) Multienvironment model: Factor analytic (reduced rank) model (FA) 10) Two stage analysis When the breeder decides to run a trial and apply selection in a single environment (whether because the amount of seed is a limitation or there's no availability for a location) the breeder takes the risk of selecting material for a target population of environments (TPEs) using an environment that is not representative of the larger TPE. Therefore, many breeding programs try to base their selection decision on multi-environment trial (MET) data. Models could be adjusted by adding additional information like spatial information, experimental design information, etc. In this tutorial we will focus mainly on the covariance structures for GxE and the incorporation of relationship matrices for the genotype effect. ## 1) Single environment model A single-environment model is the one that is fitted when the breeding program can only afford one location, leaving out the possible information available from other environments. This will be used to further expand to GxE models. ```{r} library(sommer) data(DT_example, package="enhancer") DT <- DT_example A <- A_example Ai <- solve(A) Ai <- as(as(as( Ai, "dMatrix"), "generalMatrix"), "CsparseMatrix") attr(Ai, "inverse")=TRUE ansSingle <- mmes(Yield~1, random= ~ vsm(ism(Name), Gu=Ai), rcov= ~ units, data=DT, verbose = FALSE) summary(ansSingle) ``` In this model, the random term is the germplasm effect (here called `Name`). For the sake of example, a relationship structure among the levels of `Name` is supplied through `Gu`. In the code above `Ai` is the sparse precision matrix corresponding to the relationship matrix `A`, and `attr(Ai, "inverse")=TRUE` tells sommer that the supplied matrix is already in precision form. More generally, a non-diagonal relationship structure can be used to model covariance among genotype levels. ## 2) MET: main effect model A multi-environment model is the one that is fitted when the breeding program can afford more than one location. The main effect model assumes that GxE doesn't exist and that the main genotype effect plus the fixed effect for environment is enough to predict the genotype effect in all locations of interest. ```{r} ansMain <- mmes(Yield~Env, random= ~ vsm(ism(Name), Gu=Ai), rcov= ~ units, data=DT, verbose = FALSE) summary(ansMain) ``` ## 3) MET: diagonal model (DG) A multi-environment model is fitted when observations are available from more than one environment. The diagonal GxE model allows the genetic variance to differ among environments while constraining the genetic covariance between different environments to zero. In the current `vsm()` parameterization, `dsm(Env)` is a normalized diagonal covariance factor: one product-level variance scale is estimated by `vsm()` and the remaining diagonal parameters are relative variance ratios. Thus the model has an environment-specific genetic variance without estimating cross-environment genetic covariances. The fixed environment effect and the environment-specific genotype BLUP are then combined to predict performance in each environment. ```{r} ansDG <- mmes(Yield~Env, random= ~ vsm(dsm(Env),ism(Name), Gu=Ai), rcov= ~ units, data=DT, verbose = FALSE) summary(ansDG) ``` ## 4) MET: compund symmetry model (CS) A multi-environment model is fitted when observations are available from more than one environment. In the code below, compound-symmetry-like GxE behavior is represented by two independent random terms: a genotype main effect shared across environments and an environment-by-genotype deviation. Their covariance contribution is the sum of the two random-effect covariance matrices. Under the identity structures used here this implies a common covariance between observations on the same genotype in different environments, while the GxE deviation adds environment-specific variance. The fixed environment effect, genotype main-effect BLUP, and genotype-by-environment BLUP together describe performance in each environment. ```{r} E <- diag(length(unique(DT$Env)));rownames(E) <- colnames(E) <- unique(DT$Env) Ei <- solve(E) Ai <- solve(A) EAi <- kronecker(Ei,Ai, make.dimnames = TRUE) Ei <- as(as(as( Ei, "dMatrix"), "generalMatrix"), "CsparseMatrix") Ai <- as(as(as( Ai, "dMatrix"), "generalMatrix"), "CsparseMatrix") EAi <- as(as(as( EAi, "dMatrix"), "generalMatrix"), "CsparseMatrix") attr(Ai, "inverse")=TRUE attr(EAi, "inverse")=TRUE ansCS <- mmes(Yield~Env, random= ~ vsm(ism(Name), Gu=Ai) + vsm(ism(Env:Name), Gu=EAi), rcov= ~ units, data=DT, verbose = FALSE) summary(ansCS) ``` ## 5) MET: unstructured model (US) The unstructured GxE model allows a distinct genetic variance for each environment and a distinct genetic covariance for every pair of environments. In the current `vsm()` architecture, `usm(Env)` is represented by a normalized Cholesky covariance factor and `vsm()` supplies the single product-level variance scale. This parameterization guarantees positive definiteness for admissible working parameters while retaining the flexibility of an unstructured covariance matrix. Because the number of covariance-shape parameters grows rapidly with the number of environments, these models can still be statistically weakly identified or computationally demanding. The fixed environment effect and the correlated environment-specific genotype BLUPs are used to predict performance in each environment. ```{r} ansUS <- mmes(Yield~Env, random= ~ vsm(usm(Env),ism(Name), Gu=Ai), rcov= ~ units, data=DT, verbose = FALSE) summary(ansUS) ``` ## 6) MET: random regression model A random regression model represents environmental response with a continuous covariate or a set of basis functions, here Legendre polynomials of the numeric environment index. Genotype-specific coefficients on these basis functions describe reaction norms across environments. The covariance structure assigned to the basis coefficients determines how the random intercept, slopes, and higher-order coefficients vary and covary. Consequently, the number of covariance parameters depends both on the polynomial order and on the covariance structure placed on the basis coefficients. ```{r} library(orthopolynom) DT$EnvN <- as.numeric(as.factor(DT$Env)) ansRR <- mmes(Yield~Env, random= ~ vsm(dsm(leg(EnvN,1)),ism(Name)), rcov= ~ units, data=DT, verbose = FALSE) summary(ansRR) ``` In addition, an unstructured, diagonal or other variance-covariance structure can be put on top of the polynomial model: ```{r} library(orthopolynom) DT$EnvN <- as.numeric(as.factor(DT$Env)) ansRR <- mmes(Yield~Env, random= ~ vsm(usm(leg(EnvN,1)),ism(Name)), rcov= ~ units, data=DT, verbose = FALSE) summary(ansRR) ``` ## 7) Other GxE covariance structures Many structured covariance models can be used for GxE when their assumptions are appropriate for the ordering or relationship among environments. In the example below `csm(Env)` specifies a homogeneous correlation structure for the environment covariance factor. Other available structures, such as autoregressive models for meaningfully ordered environments, can be substituted within the same `vsm()` product architecture. ```{r} ansAR1 <- mmes(Yield~Env, random= ~ vsm(csm(Env),ism(Name)), rcov= ~ units, data=DT, verbose = FALSE) summary(ansAR1) ``` ## 8) Finlay-Wilkinson regression ```{r} data(DT_h2, package="enhancer") DT <- DT_h2 ## build the environmental index ei <- aggregate(y~Env, data=DT,FUN=mean) colnames(ei)[2] <- "envIndex" ei$envIndex <- ei$envIndex - mean(ei$envIndex,na.rm=TRUE) # center the envIndex to have clean VCs ei <- ei[with(ei, order(envIndex)), ] ## add the environmental index to the original dataset DT2 <- merge(DT,ei, by="Env") # numeric by factor variables like envIndex:Name can't be used in the random part like this # they need to come with the vsm() structure DT2 <- DT2[with(DT2, order(Name)), ] mix2 <- mmes(y~ envIndex, random=~ Name + vsm(ism(envIndex),ism(Name)), data=DT2, rcov=~vsm(dsm(Name),ism(units)), nIters = 50, verbose = FALSE ) # summary(mix2)$varcomp b=mix2$uList$`vsm(ism(envIndex), ism(Name` # adaptability (b) or genotype slopes mu=mix2$uList$`vsm(ism(Name`# general adaptation (mu) or main effect e=sqrt(summary(mix2)$varcomp[-c(1:2),"estimate"]) # error variance for each individual ## general adaptation (main effect) vs adaptability (response to better environments) plot(mu[,1]~b[,1], ylab="general adaptation", xlab="adaptability") text(y=mu[,1],x=b[,1], labels = rownames(mu), cex=0.5, pos = 1) ## prediction across environments Dt <- mix2$Dtable Dt[1,"average"]=TRUE Dt[2,"include"]=TRUE Dt[3,"include"]=TRUE mix2 <- postPEV(mix2, mode=2) pp <- predict(mix2,Dtable = Dt, D="Name") preds <- pp$pvals # preds[with(preds, order(-predicted.value)), ] ## performance vs stability (deviation from regression line) plot(preds[,2]~e, ylab="performance", xlab="stability") text(y=preds[,2],x=e, labels = rownames(mu), cex=0.5, pos = 1) ``` ## 9) Factor analytic and reduced rank models When the number of environments is large, a fully unstructured genetic covariance requires \(q(q+1)/2\) covariance parameters for \(q\) environments and can become weakly identified relative to the available information. Reduced-rank and factor-analytic representations provide a more parsimonious alternative. In the first implementation below, we show the use of the factor analytic structure `fam()` where loadings and specific variances are estimated through REML and produce predictions. ```{r} data(DT_h2, package="enhancer") DT <- DT_h2 DT=DT[with(DT, order(Env)), ] head(DT) indNames <- na.omit(unique(DT$Name)) A <- diag(length(indNames)) rownames(A) <- colnames(A) <- indNames # factor analytic model with 2 factors ansFA2 <- mmes(y~Env, random=~vsm( fam(Env, 2) , ism(Name)) , rcov=~units, nIters = 100, verbose = FALSE, data=DT) Dt <- ansFA2$Dtable; Dt Dt[1:2,"average"]=TRUE Dt[3,c("average","include")]=TRUE ppfa <- predict(ansFA2, D="Env:Name", Dtable = Dt) head(ppfa$pvals) ``` If needed, loadings and factor scores can be extracted with the `loadings_mmes()` and `scores_mmes()` functions. The second implementation below uses the reduced rank implementation available with the `rrm()` function to estimate loadings through REML and we use a diagonal term (dsm) to fit the environment specific variances. ```{r} # reduced rank model with 2 factors ansRR2 <- mmes(y~Env, henderson=TRUE, random=~vsm( rrm(Env, 2) , ism(Name)) + # rr vsm(dsm(Env), ism(Name)), # diag rcov=~units, nIters = 100, verbose = FALSE, data=DT) Dt <- ansRR2$Dtable; Dt Dt[1:2,"average"]=TRUE Dt[3:4,c("average","include")]=TRUE pprr <- predict(ansRR2, D="Env:Name", Dtable = Dt) head(pprr$pvals) # compare plot(ppfa$pvals[,"predicted.value"],pprr$pvals[,"predicted.value"]) ``` Last but not least, it is possible to avoid the estimation of loadings by defining fixed loadings as covariates of a random regression coming from a diagonal model that uses the GxE predictions to derive Sigma and the fixed loadings. ```{r} # fit diagonal model first to produce H matrix ansDG <- mmes(y~Env, henderson=TRUE, random=~ vsm(dsm(Env), ism(Name)), rcov=~units, nIters = 100, data=DT, verbose = FALSE) H0 <- ansDG$uList$`vsm(dsm(Env), ism(Name))` # GxE table # # reduced rank model # ansFA <- mmes(y~Env, henderson=TRUE, # random=~vsm( usm(rrmat(Env, H = H0, nPC = 2)) , ism(Name)) + # rr # vsm(dsm(Env), ism(Name)), # diag # rcov=~units, # # we recommend giving more iterations to these models # nIters = 100, verbose = FALSE, # # we recommend giving more EM iterations at the beggining # data=DT) # # vcFA <- ansFA$theta[[1]] # vcDG <- ansFA$theta[[2]] # # loadings=with(DT, rrmat(Env, nPC = 2, H = H0, returnGamma = TRUE) )$Gamma # scores <- ansFA$uList[[1]] # # vcUS <- loadings %*% vcFA %*% t(loadings) # G <- vcUS + vcDG # # colfunc <- colorRampPalette(c("steelblue4","springgreen","yellow")) # # hv <- heatmap(cov2cor(G), col = colfunc(100), symm = TRUE) # # uFA <- scores %*% t(loadings) # uDG <- ansFA$uList[[2]] # u <- uFA + uDG # # plot(ppfa$pvals[,"predicted.value"],as.vector(t(u))) # plot(as.vector(t(u)),pprr$pvals[,"predicted.value"]) ``` ## 10) Two stage analysis In two-stage analyses, a first-stage mixed model is commonly fitted to account for experimental-design and field variation while treating the entry effects of interest as fixed, producing adjusted entry means (BLUEs or EMMs). Their estimation-error covariance should be carried into the second stage rather than treating the adjusted means as independent observations with equal precision. In the example below, `computeCi = 2` requests the complete coefficient-matrix inverse after fitting each first-stage model; the fixed-effect block is used to construct the corresponding precision contribution for the second-stage weight matrix `W`. The second-stage mixed model then analyzes the adjusted means while retaining this first-stage precision information. ```{r} ########## ## stage 1 ## use mmes for dense field trials ########## data(DT_h2, package="enhancer") DT <- DT_h2 head(DT) envs <- unique(DT$Env) BLUEL <- list() XtXL <- list() for(i in 1:length(envs)){ ans1 <- mmes(y~Name-1, random=~Block, verbose=FALSE, computeCi = 2, data=droplevels(DT[which(DT$Env == envs[i]),]) ) ans1$Beta$Env <- envs[i] BLUEL[[i]] <- data.frame( Effect=factor(rownames(ans1$b)), Estimate=ans1$b[,1], Env=factor(envs[i])) # to be comparable to 1/(se^2) = 1/PEV = 1/Ci = 1/[(X'X)inv] XtXL[[i]] <- solve(ans1$Ci[1:nrow(ans1$b),1:nrow(ans1$b)]) } DT2 <- do.call(rbind, BLUEL) OM <- Reduce(adiag1,lapply(XtXL,as.matrix)) ########## ## stage 2 ## use mmes for sparse equation ########## m <- matrix(1/var(DT2$Estimate, na.rm = TRUE)) ans2 <- mmes(Estimate~Env, henderson=TRUE, random=~ Effect + Env:Effect, rcov = ~ vsm( ism(units), sigma2 = 1, fixedSigma2 = TRUE ), W=OM, verbose=FALSE, data=DT2 ) summary(ans2)$varcomp ``` ## Literature Covarrubias-Pazaran G. 2016. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6):1-15. Covarrubias-Pazaran G. 2018. Software update: Moving the R package sommer to multivariate mixed models for genome-assisted prediction. doi: https://doi.org/10.1101/354639 Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp. Gilmour et al. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51(4):1440-1450. Henderson C.R. 1975. Best Linear Unbiased Estimation and Prediction under a Selection Model. Biometrics vol. 31(2):423-447. Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723. Lee, D.-J., Durban, M., and Eilers, P.H.C. (2013). Efficient two-dimensional smoothing with P-spline ANOVA mixed models and nested bases. Computational Statistics and Data Analysis, 61, 22 - 37. Lee et al. 2015. MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Cold Spring Harbor. doi: http://dx.doi.org/10.1101/027201. Maier et al. 2015. Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder. Am J Hum Genet; 96(2):283-294. Rodriguez-Alvarez, Maria Xose, et al. Correcting for spatial heterogeneity in plant breeding experiments with P-splines. Spatial Statistics 23 (2018): 52-71. Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco. Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208. Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.