reduced rank matrix
rrm creates a reduced rank factor analytic matrix by selecting the n vectors of the L matrix of the Cholesky decomposition or the U vectors of the SVD decomposition (loadings or latent covariates) to create a new incidence matrix of latent covariates that can be used with the lmebreed solver to fit random regressions on the latent covariates.
rrm(x=NULL, H=NULL, nPC=2, returnGamma=FALSE, cholD=TRUE)
x: vector of the dataset containing the variable to be used to form the incidence matrix.H: two-way table of identifiers (rows; e.g., genotypes) by features (columns; e.g., environments) effects. Row names and column names are required. No missing data is allowed.nPC: number of principal components to keep from the loadings matrix.returnGamma: a TRUE/FALSE argument specifying if the function should return the matrix of loadings used to build the incidence matrix for the model. The default is FALSE so it returns only the incidence matrix.cholD: a TRUE/FALSE argument specifying if the Cholesky decomposition should be calculated or the singular value decomposition should be used instead.This implementation of a version of the reduced rank factor analytic models uses the so-called principal component (PC) models (Meyer, 2009) which assumes specific effects (psi) are equal to 0. The model is as follows:
y = Xb + Zu + e
where the variance of u ~ MVN(0, Sigma)
Sigma = (Gamma_t Gamma) + Psi
Extended factor analytic model :
y = Xb + Z(I Gamma)c + Zs + e = Xb + Z*c + Zs + e
where y is the response variable, X and Z are incidence matrices for fixed and random effects respectively, I is a diagonal matrix, Gamma are the factor loadings for c common factor scores, and s are the specific effects, e is the vector of residuals.
Reduced rank model:
y = Xb + Z(I Gamma)c + e = Xb + Z*c + e
which is equal to the one above but assumes specific effects = 0.
The algorithm in rrm is the following:
uses a wide-format table of timevar (m columns) by idvar (q rows) named H to form the initial variance-covariance matrix (Sigma) which is calculated as Sigma = H'H of dimensions m x m (column dimensions, e.g., environments x environments).
The Sigma matrix is then center and scaled.
A Cholesky (L matrix) or SVD decomposition (U D V') is performed in the Sigma matrix.
n vectors from L (when Cholesky is used) or U sqrt(D) (when SVD is used) are kept to form Gamma. Gamma = L[,1:nPc] or Gamma = U[,1:nPC]. These are the so-called loadings (L for all loadings, Gamma for the subset of loadings).
Gamma is used to form a new incidence matrix as Z* = Z Gamma
This matrix is later used for the REML machinery to be used with the usc (unstructured) or smm (diagonal) structures to estimate variance components and factor scores. The resulting BLUPs from the mixed model are the optimized factor scores. Pretty much as a random regression over latent covariates.
This implementation does not update the loadings (latent covariates) during the REML process, only estimates the REML factor scores for fixed loadings. This is different to other software (e.g., asreml) where the loadings are updated during the REML process as well.
BLUPs for genotypes in all locations can be recovered as:
u = Gamma * u_scores
The resulting loadings (Gamma) and factor scores can be thought as an equivalent to the classical factor analysis.
As an additional information, notice that we calculate the factor loadings from BLUPs and the mixed model only calculates the factor scores. This is different to the asreml software where loadings are calculated as variance components through REML. Despite the difference we have run multiple examples and simulations and the BLUPs from both approaches are on average >0.98 correlated so you can be confident that our approach is robust.
Giovanny Covarrubias-Pazaran (2024). lme4breeding: enabling genetic evaluation in the age of genomic data. To be submitted to Bioinformatics.
Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48.
Meyer K (2009) Factor analytic models for genotype by environment type problems and structured covariance matrices. Genetics Selection Evolution, 41:21
Giovanny Covarrubias-Pazaran
data(DT_h2) DT <- DT_h2 DT=DT[with(DT, order(Env)), ] indNames <- na.omit(unique(DT$Name)) A <- diag(length(indNames)) rownames(A) <- colnames(A) <- indNames # fit diagonal model first to produce H matrix Z <- with(DT, smm(Env)) diagFormula <- paste0( "y ~ Env + (0+", paste(colnames(Z), collapse = "+"), "|| Name)") for(i in 1:ncol(Z)){DT[,colnames(Z)[i]] <- Z[,i]} print(as.formula(diagFormula)) ans1a <- lmebreed(as.formula(diagFormula), relmat = list(Name=A), data=DT) vc <- VarCorr(ans1a); print(vc,comp=c("Variance")) H0 <- ranef(ans1a)$Name # GxE table # reduced rank model Z <- with(DT, rrm(Env, H = H0, nPC = 3)) Zd <- with(DT, smm(Env)) faFormula <- paste0( "y ~ Env + (0+", paste(colnames(Z), collapse = "+"), "| Name) + (0+",paste(colnames(Zd), collapse = "+"), "|| Name)") for(i in 1:ncol(Z)){DT[,colnames(Z)[i]] <- Z[,i]} print(as.formula(faFormula)) ansFA <- lmebreed(as.formula(faFormula), relmat = list(Name=A), data=DT) vc <- VarCorr(ansFA); print(vc,comp=c("Variance")) ve <- attr(vc, "sc")^2; ve loadings=with(DT, rrm(Env, nPC = 3, H = H0, returnGamma = TRUE) )$Gamma Gint <- loadings %*% vc$Name %*% t(loadings) Gspec <- diag( unlist(lapply(vc[2:16], function(x){x[[1]]})) ) G <- Gint + Gspec # lattice::levelplot(cov2cor(G)) # colfunc <- colorRampPalette(c("steelblue4","springgreen","yellow")) # hv <- heatmap(cov2cor(G), col = colfunc(100), symm = TRUE) u <- ranef(ansFA)$Name uInter <- as.matrix(u[,1:3]) %*% t(as.matrix(loadings)) uSpec <- as.matrix(u[,-c(1:3)]) u <- uSpec + uInter
The lmebreed solver.
Useful links