Fitting genotype by environment models in sommer

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:

Single environment model
Multienvironment model: Main effect model
Multienvironment model: Diagonal model (DG)
Multienvironment model: Compund symmetry model (CS)
Multienvironment model: Unstructured model (US)
Multienvironment model: Random regression model (RR)
Multienvironment model: Other covariance structures for GxE
Multienvironment model: Finlay-Wilkinson regression
Multienvironment model: Factor analytic (reduced rank) model (FA)
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.

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,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansSingle)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -85.84241 173.6848 176.9052     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                      term parameter estimate StdError Zratio
## 1 vsm(ism(Name), Gu = Ai)    sigma2    10.22    3.219  3.174
## 2         vsm(ism(units))    sigma2    10.22    1.031  9.912
## ============================================================
## Fixed effects:
##           Estimate Std.Error t.value
## Intercept    11.73        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

ansMain <- mmes(Yield~Env,
  random= ~ vsm(ism(Name), Gu=Ai),
  rcov= ~ units,
  data=DT,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansMain)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -40.69566 87.39131 97.05238     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                      term parameter estimate StdError Zratio
## 1 vsm(ism(Name), Gu = Ai)    sigma2    6.631   2.0222  3.279
## 2         vsm(ism(units))    sigma2    6.631   0.7086  9.358
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.337        NA      NA
## EnvCA.2012   -5.625        NA      NA
## EnvCA.2013   -6.140        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

ansDG <- mmes(Yield~Env,
  random= ~ vsm(dsm(Env),ism(Name), Gu=Ai),
  rcov= ~ units,
  data=DT,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansDG)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -31.75412 69.50824 79.16931     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                                term         parameter estimate StdError Zratio
## 1 vsm(dsm(Env), ism(Name), Gu = Ai) variance[CA.2011]   10.058    3.467  2.901
## 2 vsm(dsm(Env), ism(Name), Gu = Ai) variance[CA.2012]    4.649    2.037  2.282
## 3 vsm(dsm(Env), ism(Name), Gu = Ai) variance[CA.2013]    5.187    2.054  2.525
## 4                   vsm(ism(units))            sigma2    6.631    1.211  5.476
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.628        NA      NA
## EnvCA.2012   -5.979        NA      NA
## EnvCA.2013   -6.676        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

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,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))

## Adding 29 additional Gu levels to the main-effect model matrix: CA.2013:CO05061-2P, CA.2013:AC01151-5W, CA.2013:MSS165-2Y, CA.2013:NY148, CA.2013:W8539-2Y, CA.2013:W8615-5, CA.2013:MSR148-4, CA.2013:CO00270-7W, CA.2011:Manistee(MSL292-A), CA.2011:AC05153-1W ...

summary(ansCS)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -29.94453 65.88906 75.55013     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                           term parameter estimate StdError Zratio
## 1      vsm(ism(Name), Gu = Ai)    sigma2    3.316    1.396  2.375
## 2 vsm(ism(Env:Name), Gu = EAi)    sigma2    3.316    1.392  2.382
## 3              vsm(ism(units))    sigma2    6.631    1.211  5.476
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.490        NA      NA
## EnvCA.2012   -5.789        NA      NA
## EnvCA.2013   -6.380        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

ansUS <- mmes(Yield~Env,
  random= ~ vsm(usm(Env),ism(Name), Gu=Ai),
  rcov= ~ units,
  data=DT,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansUS)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -30.01355 66.02711 75.68817     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                                term                   parameter estimate
## 1 vsm(usm(Env), ism(Name), Gu = Ai)           variance[CA.2011]  10.0577
## 2 vsm(usm(Env), ism(Name), Gu = Ai) covariance[CA.2012,CA.2011]   0.6838
## 3 vsm(usm(Env), ism(Name), Gu = Ai) covariance[CA.2013,CA.2011]   0.7223
## 4 vsm(usm(Env), ism(Name), Gu = Ai)           variance[CA.2012]   4.6492
## 5 vsm(usm(Env), ism(Name), Gu = Ai) covariance[CA.2013,CA.2012]   0.4911
## 6 vsm(usm(Env), ism(Name), Gu = Ai)           variance[CA.2013]   5.1873
## 7                   vsm(ism(units))                      sigma2   6.6314
##   StdError Zratio
## 1    3.900 2.5791
## 2    2.359 0.2899
## 3    2.468 0.2927
## 4    2.143 2.1692
## 5    1.596 0.3078
## 6    2.152 2.4104
## 7    1.211 5.4764
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.595        NA      NA
## EnvCA.2012   -5.936        NA      NA
## EnvCA.2013   -6.609        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

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,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansRR)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -42.54287 91.08575 100.7468     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                                term      parameter estimate StdError Zratio
## 1 vsm(dsm(leg(EnvN, 1)), ism(Name)) variance[leg0]    5.236   1.5790  3.316
## 2 vsm(dsm(leg(EnvN, 1)), ism(Name)) variance[leg1]    8.026   3.7440  2.144
## 3                   vsm(ism(units))         sigma2    6.631   0.9357  7.087
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.617        NA      NA
## EnvCA.2012   -5.936        NA      NA
## EnvCA.2013   -6.704        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

In addition, an unstructured, diagonal or other variance-covariance structure can be put on top of the polynomial model:

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,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansRR)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -43.42522 92.85043 102.5115     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                                term             parameter estimate StdError
## 1 vsm(usm(leg(EnvN, 1)), ism(Name))        variance[leg0]   5.2364   1.6020
## 2 vsm(usm(leg(EnvN, 1)), ism(Name)) covariance[leg1,leg0]   0.6483   1.6145
## 3 vsm(usm(leg(EnvN, 1)), ism(Name))        variance[leg1]   8.0264   3.8696
## 4                   vsm(ism(units))                sigma2   6.6314   0.9337
##   Zratio
## 1 3.2686
## 2 0.4016
## 3 2.0742
## 4 7.1021
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.627        NA      NA
## EnvCA.2012   -5.945        NA      NA
## EnvCA.2013   -6.713        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

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.

ansAR1 <- mmes(Yield~Env,
  random= ~ vsm(csm(Env),ism(Name)),
  rcov= ~ units,
  data=DT,
  control = mmes_control(convergence = list(maxIters = 1L), diagnostics = list(verbose = FALSE)))
summary(ansAR1)

## ============================================================
##          Multivariate Linear Mixed Model fit by  REML         
## **********************  sommer 4.4  ********************** 
## ============================================================
##          logLik      AIC      BIC Method Converge
## Value -32.62932 71.25865 80.91972     AI    FALSE
## ============================================================
## Variance-Covariance components:
##                       term parameter estimate StdError Zratio
## 1 vsm(csm(Env), ism(Name))  variance    6.631   1.6116 4.1149
## 2 vsm(csm(Env), ism(Name))       rho    0.100   0.1705 0.5865
## 3          vsm(ism(units))    sigma2    6.631   1.2123 5.4701
## ============================================================
## Fixed effects:
##            Estimate Std.Error t.value
## Intercept    16.605        NA      NA
## EnvCA.2012   -5.936        NA      NA
## EnvCA.2013   -6.612        NA      NA
## ============================================================
## Use the '$' sign to access results and parameters

8) Finlay-Wilkinson regression

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)),
  control = mmes_control(
    # convergence = list(maxIters = 50),
  diagnostics = list(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.

data(DT_h2, package="enhancer")
DT <- DT_h2
DT=DT[with(DT, order(Env)), ]
head(DT)

##          Name     Env Loc Year     Block  y
## 67   MSL007-B CA.2011  CA 2011 CA.2011.2  5
## 105  MSL007-B CA.2011  CA 2011 CA.2011.1  6
## 308  MSK061-4 CA.2011  CA 2011 CA.2011.2  9
## 393  MSK061-4 CA.2011  CA 2011 CA.2011.1 10
## 469 MSR169-8Y CA.2011  CA 2011 CA.2011.1 11
## 471     NY148 CA.2011  CA 2011 CA.2011.1 11

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,
  data=DT,
  control = mmes_control(
    convergence = list(maxIters = 1L),
  diagnostics = list(verbose = FALSE)))

Dt <- ansFA2$Dtable; Dt

##     type                        term include average levels
## 1  fixed                           1   FALSE   FALSE   NULL
## 2  fixed                         Env   FALSE   FALSE   NULL
## 3 random vsm(fam(Env, 2), ism(Name))   FALSE   FALSE   NULL

Dt[1:2,"average"]=TRUE
Dt[3,c("average","include")]=TRUE

ppfa <- predict(ansFA2, D="Env:Name", Dtable = Dt)
head(ppfa$pvals)

##                            Env:Name predicted.value std.error      pev
## CA.2011:A01143-3C CA.2011:A01143-3C        16.24857  1.339181 1.793405
## CA.2012:A01143-3C CA.2012:A01143-3C        11.68609  1.318696 1.738960
## CA.2013:A01143-3C CA.2013:A01143-3C        17.45958  1.196669 1.432017
## FL.2011:A01143-3C FL.2011:A01143-3C        11.73320  1.204849 1.451660
## FL.2012:A01143-3C FL.2012:A01143-3C        12.10790  1.194296 1.426343
## FL.2013:A01143-3C FL.2013:A01143-3C        12.36608  1.742795 3.037336
##                        var.u
## CA.2011:A01143-3C 3.26201872
## CA.2012:A01143-3C 3.38728599
## CA.2013:A01143-3C 1.64188386
## FL.2011:A01143-3C 1.61499095
## FL.2012:A01143-3C 1.64921738
## FL.2013:A01143-3C 0.04745611

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.

# reduced rank model with 2 factors
ansRR2 <- mmes(y~Env,
  random=~vsm( rrm(Env, 2) , ism(Name)) + # rr
                vsm(dsm(Env), ism(Name)),
  # diag
              rcov=~units,
  data=DT,
  control = mmes_control(formulation = "henderson",
              convergence = list(maxIters = 1L),
  diagnostics = list(verbose = FALSE)))


Dt <- ansRR2$Dtable; Dt

##     type                        term include average levels
## 1  fixed                           1   FALSE   FALSE   NULL
## 2  fixed                         Env   FALSE   FALSE   NULL
## 3 random vsm(rrm(Env, 2), ism(Name))   FALSE   FALSE   NULL
## 4 random    vsm(dsm(Env), ism(Name))   FALSE   FALSE   NULL

Dt[1:2,"average"]=TRUE
Dt[3:4,c("average","include")]=TRUE

pprr <- predict(ansRR2, D="Env:Name", Dtable = Dt)
head(pprr$pvals)

##                            Env:Name predicted.value std.error       pev
## CA.2011:A01143-3C CA.2011:A01143-3C        16.66591 1.4532657 2.1119813
## CA.2012:A01143-3C CA.2012:A01143-3C        11.70699 1.3041714 1.7008630
## CA.2013:A01143-3C CA.2013:A01143-3C        17.22872 1.1724072 1.3745386
## FL.2011:A01143-3C FL.2011:A01143-3C        11.94660 0.9494034 0.9013668
## FL.2012:A01143-3C FL.2012:A01143-3C        12.21767 0.7977691 0.6364355
## FL.2013:A01143-3C FL.2013:A01143-3C        12.31767 0.9682978 0.9376006
##                        var.u
## CA.2011:A01143-3C 5.64151020
## CA.2012:A01143-3C 3.32199017
## CA.2013:A01143-3C 1.47782119
## FL.2011:A01143-3C 0.50308062
## FL.2012:A01143-3C 0.25959742
## FL.2013:A01143-3C 0.06361307

# 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.

# fit diagonal model first to produce H matrix
# ansDG <- mmes(y~Env,
#   random=~ vsm(dsm(Env), ism(Name)),
#   rcov=~units,
#   data=DT,
#   control = mmes_control(formulation = "henderson",
#   # convergence = list(maxIters = 100),
#   diagnostics = list(verbose = FALSE)))
# 
# H0 <- ansDG$uList$`vsm(dsm(Env), ism(Name))` # GxE table

# # reduced rank model
# ansFA <- mmes(y~Env,
#               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
#               control=mmes_control(formulation="henderson",
#                 convergence=list(maxIters=100),
#                 diagnostics=list(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, control$output$inverse = 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.

##########
## stage 1
## use mmes for dense field trials
##########
data(DT_h2, package="enhancer")
DT <- DT_h2
head(DT)

##                 Name     Env Loc Year     Block y
## 1            W8822-3 FL.2012  FL 2012 FL.2012.1 2
## 2            W8867-7 FL.2012  FL 2012 FL.2012.2 2
## 3           MSL007-B MO.2011  MO 2011 MO.2011.1 3
## 4         CO00270-7W FL.2012  FL 2012 FL.2012.2 3
## 5 Manistee(MSL292-A) FL.2013  FL 2013 FL.2013.2 3
## 6           MSM246-B FL.2012  FL 2012 FL.2012.2 3

envs <- unique(DT$Env)
BLUEL <- list()
XtXL <- list()
for(i in 1:length(envs)){
  ans1 <- mmes(y~Name-1,
  random=~Block,
  data=droplevels(DT[which(DT$Env == envs[i]),]),
  control = mmes_control(convergence = list(maxIters = 1L),
  diagnostics = list(verbose = FALSE),
  output = list(inverse = 2)))
  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,
  random=~ Effect + Env:Effect,
  rcov = ~ vsm(
               ism(units),
               sigma2 = 1,
               fixedSigma2 = TRUE
             ),
  W=OM,
  data=DT2,
  control = mmes_control(formulation = "henderson",
  convergence = list(maxIters = 1L),
  diagnostics = list(verbose = FALSE)))
summary(ans2)$varcomp

##                                              term parameter estimate StdError
## 1                                vsm(ism(Effect))    sigma2 2.432911 0.684542
## 2                            vsm(ism(Env:Effect))    sigma2 4.262498 0.384575
## 3 vsm(ism(units), sigma2 = 1, fixedSigma2 = TRUE)    sigma2 1.000000 0.000000
##      Zratio
## 1  3.554071
## 2 11.083657
## 3        NA

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.