R/H_matrix.R

In ASRgenomics: Complementary Genomic Functions

#' Generates the hybrid \eqn{H} matrix
#'
#' The single-step GBLUP approach combines the information from the pedigree relationship matrix
#' \eqn{\boldsymbol{A}} and the genomic relationship matrix \eqn{\boldsymbol{G}} in one
#' hybrid relationship matrix called \eqn{\boldsymbol{H}}.
#' This function will calculate directly this matrix \eqn{\boldsymbol{H}}.
#' The user should provide the matrices \eqn{\boldsymbol{A}} or
#' its inverse (only one of these is required) and the
#' inverse of the matrix \eqn{\boldsymbol{G}} (\eqn{\boldsymbol{G_{inv}}}) in its full form.
#' Individual names should
#' be assigned to  \code{rownames} and \code{colnames}, and individuals from
#' \eqn{\boldsymbol{G_{inv}}} are verified to be all a subset within individuals from
#' \eqn{\boldsymbol{A}} (or \eqn{\boldsymbol{A_{inv}}}).
#' This function is a wrapper of  the \link{H.inverse} function.
#'
#' @param A Input of the pedigree relationship matrix \eqn{\boldsymbol{A}}
#' in full form (\eqn{na \times na}) (default = \code{NULL}).
#' @param Ainv Input of the inverse of the pedigree relationship matrix
#' \eqn{\boldsymbol{A}^{-1}} in full form (\eqn{na \times na}) (default = \code{NULL}).
#' @param Ginv Input of the inverse of the genomic relationship matrix
#' \eqn{\boldsymbol{G}^{-1}} in full form (\eqn{ng \times ng}) (default = \code{NULL}).
#' @param lambda The scaling factor for \eqn{(\boldsymbol{G}^{-1}-\boldsymbol{A}^{-1}_{22})} (default = \code{NULL}).
#' @param tau The scaling factor for \eqn{\boldsymbol{G}^{-1}} (default = \code{1}).
#' @param omega The scaling factor for \eqn{\boldsymbol{A}^{-1}_{22}} (default = \code{1}).
#' @param sparseform If \code{TRUE} it generates the requested matrix in sparse form to be used
#' directly in \pkg{asreml} with required attributes (default = \code{FALSE}).
#' @param keep.order If \code{TRUE} the original order of the individuals from the
#' \eqn{\boldsymbol{A}} or \eqn{\boldsymbol{A_{inv}}} matrix is kept.
#' Otherwise the non-genotyped individuals are placed first and
#' then genotyped individuals (default = \code{TRUE}).
#' @param digits Set up the number of digits used to round the output matrix (default = \code{8}).
#' @param message If \code{TRUE} diagnostic messages are printed on screen (default = \code{TRUE}).
#'
#' @return
#' The hybrid matrix \eqn{\boldsymbol{H}} matrix, in full or sparse form.
#'
#' @md
#' @details
#' This function is currently equivalent to using \link{H.inverse} with (\code{inverse = FALSE}).
#'
#' \if{html}{
#' The \eqn{\boldsymbol{H}} matrix is obtained with the following equations:
#'   \deqn{\boldsymbol{H}=\boldsymbol{A}+\begin{bmatrix}\boldsymbol{A}_{12}\boldsymbol{A}_{22}^{-1}(\boldsymbol{G}-\boldsymbol{A}_{22})\boldsymbol{A}_{22}^{-1}\boldsymbol{A}_{21}&\boldsymbol{A}_{12}\boldsymbol{A}_{22}^{-1}(\boldsymbol{G}-\boldsymbol{A}_{22})\\(\boldsymbol{G}-\boldsymbol{A}_{22})\boldsymbol{A}_{22}^{-1}\boldsymbol{A}_{21}&(\boldsymbol{G}-\boldsymbol{A}_{22})\end{bmatrix}}
#' }
#'
#'
#' @references
#' Christensen, O.F., Lund, M.S. 2010. Genomic prediction matrix when some animals
#' are not genotyped. Gen. Sel. Evol. 42(2):1–8.
#'
#' Christensen, O., Madsen, P., Nielsen, B., Ostersen, T., and Su, G. 2012. Single-step methods
#' for genomic evaluation in pigs. Animal 6(10):1565–1571.
#'
#' Legarra, A., Aguilar, I., and Misztal, I. 2009. A relationship matrix including full
#' pedigree and genomic information. J. Dairy Sci. 92:4656-4663.
#'
#' Martini, J.W.R., Schrauf, M.F., Garcia-Baccino, C.A., Pimentel, E.C.G., Munilla, S.,
#' Rogberg-Muñoz, A., Cantet, R.J.C., Reimer, C., Gao, N., Wimmer, V., and Simianer, H. 2018.
#' The effect of the \eqn{H^{-1}} scaling factors \eqn{\tau} and \eqn{\omega}
#' on the structure of \eqn{H} in the single-step procedure.
#' Genet. Sel. Evol. 50:1-9.
#'
#' @export
#'
#' @examples
#' \donttest{
#' # Get A matrix.
#' A <- AGHmatrix::Amatrix(data = ped.pine)
#' A[1:5,1:5]
#' dim(A)
#'
#' # Read and filter genotypic data.
#' M.clean <- qc.filtering(
#'  M = geno.pine655,
#'  maf = 0.05,
#'  marker.callrate = 0.2, ind.callrate = 0.20,
#'  na.string = "-9",
#'  plots = FALSE)$M.clean
#'
#' # Get G matrix.
#' G <- G.matrix(M = M.clean, method = "VanRaden", na.string = "-9")$G
#' G[1:5, 1:5]
#' dim(G)
#'
#' # Match G2A.
#' check <- match.G2A(
#'  A = A, G = G,
#'  clean = TRUE, ord = TRUE, mism = TRUE, RMdiff = TRUE)
#'
#' # Align G matrix with A.
#' G_align <- G.tuneup(G = check$Gclean, A = check$Aclean, align = TRUE, sparseform = FALSE)$Gb
#'
#' # Get Ginverse using the aligned G.
#' Ginv <- G.inverse(G = G_align, sparseform = FALSE)$Ginv
#' Ginv[1:5, 1:5]
#' dim(Ginv)
#'
#' # Obtaining H.
#' H <- H.matrix(A = A, G = Ginv, lambda = 0.90, sparseform = FALSE)
#' H[1:5, 1:5]
#' }
#'

H.matrix <- function(A = NULL, Ainv = NULL, Ginv = NULL,
                     lambda = NULL,tau=1, omega = 1,
                     sparseform = FALSE, keep.order = TRUE,
                     digits = 8, message = TRUE){

  return(
    H.inverse(A = A, Ainv = Ainv, Ginv = Ginv,
              lambda = lambda,tau = tau, omega = omega,
              sparseform = sparseform, keep.order = keep.order,
              digits = digits, message = message,
              inverse = FALSE # This is the only one locked.
    )
  )
}