R/RcppExports.R

In statgenQTLxT: Multi-Trait and Multi-Trial Genome Wide Association Studies

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Update W and P in EMFA algorithm for homogeneous variance.
#'
#' Update W and P used in the iteration process in the EMFA algorithm in case
#' the variance is homogeneous.
#'
#' @inheritParams updateFA
#'
#' @param s A p x p sample covariance matrix.
#' @param m An integer. The order of the model.
#' @param maxDiag A numerical value for the maximum value of sigma2.
#'
#' @keywords internal
#'
updateFAHomVar <- function(s, wNew, pNew, m, maxDiag = 1e4) {
    invisible(.Call(`_statgenQTLxT_updateFAHomVar`, s, wNew, pNew, m, maxDiag))
}

#' Update W and P in EMFA algorithm
#'
#' Update W and P used in the iteration process in the EMFA algorithm.
#'
#' @param y An n x p matrix or data.frame.
#' @param wStart A p x p matrix or data.frame containing starting values for W.
#' @param m0 An integer. The order of the model.
#' @param pStart A p x p matrix or data.frame containing starting values for P.
#' @param hetVar Should an extra diagonal part be added in the model for the
#' precision matrix?
#' @param maxDiag A numerical value for the maximum value of the diagonal of P.
#' @param tolerance A numerical value. The iterating process stops if the sum
#' of the difference for P and W between two steps gets lower than this value.
#' @param maxIter A numerical value for the maximum number of iterations.
#' @param printProgress Should progress be printed during iterations?
#'
#' @keywords internal
#'
updateFA <- function(y, wStart, pStart, wNew, pNew, m0, hetVar = FALSE, maxDiag = 1e4, tolerance = 1e-4, maxIter = 100L, printProgress = FALSE) {
    invisible(.Call(`_statgenQTLxT_updateFA`, y, wStart, pStart, wNew, pNew, m0, hetVar, maxDiag, tolerance, maxIter, printProgress))
}

#' Helper function for updating precision matrix.
#'
#' Helper function for updating the precision matrices in the EMFA algorithm.
#'
#' @param m An integer, the order of the model.
#' @param nc An integer, the number of traits or genotypes.
#' @param omega A computed matrix for the current step in the algoritm.
#' @param w A model matrix for the current step in the algorithm.
#' @param p A model matrix for the current step in the algorithm.
#' @param wNew A pointer to the updated model matrix for w.
#' @param pNew A pointer to the updated model matrix for p.
#' @param cNew A pointer to the updated matrix c.
#' @param het Should an extra diagonal part be added in the model for the
#' precision matrix.
#' @param maxDiag A numerical value for the maximum value of sigma2.
#'
#' @keywords internal
#'
updatePrec <- function(m, nc, omega, w, p, wNew, pNew, cNew, het, maxDiag) {
    invisible(.Call(`_statgenQTLxT_updatePrec`, m, nc, omega, w, p, wNew, pNew, cNew, het, maxDiag))
}

#' Helper functions for the penalized EM algorithm
#'
#' \code{vecInvDiag} is a helper function for quickly computing
#' \eqn{(I + x \otimes y)^{-1}},
#' \code{tracePInvDiag} for quickly computing column sums of
#' \eqn{(I + x \otimes y)^{-1}}. Both are used in the penalized EM algorithm.
#'
#' @param x A numeric vector
#' @param y A numeric vector
#'
#' @return for \code{vecInvDiag} a matrix defined by
#' \eqn{(I + x \otimes y)^{-1}}, for \code{tracePInvDiag} a vector containing
#' the column sums of \eqn{(I + x \otimes y)^{-1}}.
#'
#' @keywords internal
#'
vecInvDiag <- function(x, y) {
    .Call(`_statgenQTLxT_vecInvDiag`, x, y)
}

#' @rdname vecInvDiag
#'
tracePInvDiag <- function(x, y) {
    .Call(`_statgenQTLxT_tracePInvDiag`, x, y)
}

#' Factor analytic variation of EM algoritm
#'
#' Implementation of the factor analytic variation of the EM algoritm as
#' proposed by Dahl et al. (2013).
#'
#' @param y An n x p matrix of observed phenotypes, on p traits or environments
#' for n individuals. No missing values are allowed.
#' @param k An n x n kinship matrix.
#' @param size_param_x An n x c covariate matrix, c being the number of
#' covariates and n being the number of genotypes. c has to be at least one
#' (typically an intercept). No missing values are allowed. If not provided a
#' vector of 1s is used.
#' @param cmHet Should an extra diagonal part be added in the model for the
#' precision matrix Cm?
#' @param dmHet Should an extra diagonal part be added in the model for the
#' precision matrix Dm?
#' @param tolerance A numerical value. The iterating process stops if the
#' difference in conditional log-likelihood between two consecutive iterations
#' drops below tolerance.
#' @param maxIter A numerical value for the maximum number of iterations.
#' @param size_param_cmStart A p x p matrix containing starting values for the
#' precision matrix Cm.
#' @param size_param_dmStart A p x p matrix containing starting values for the
#' precision matrix Dm.
#' @param mG An integer. The order of the genetic part of the model.
#' @param mE An integer. The order of the environmental part of the model.
#' @param maxDiag A numical value. The maximal value of the diagonal elements
#' in the precision matrices Cm and Dm (ignoring the low-rank part W W^t)
#' @param stopIfDecreasing Should the iterating process stop if after 50
#' iterations the log-likelihood decreases between two consecutive iterations?
#'
#' @return A list containing the following components
#' \itemize{
#' \item{\code{Vg} The genetic variance components matrix.}
#' \item{\code{Ve} The environmental variance components matrix.}
#' }
#'
#' @references Dahl et al. (2013). Network inference in matrix-variate Gaussian
#' models with non-independent noise. arXiv preprint arXiv:1312.1622.
#' @references Zhou, X. and Stephens, M. (2014). Efficient multivariate linear
#' mixed model algorithms for genome-wide association studies. Nature Methods,
#' February 2014, Vol. 11, p. 407–409
#'
#' @keywords internal
#'
EMFA <- function(y, k, size_param_x = NULL, cmHet = TRUE, dmHet = TRUE, tolerance = 1e-6, maxIter = 300L, size_param_cmStart = NULL, size_param_dmStart = NULL, mG = 1L, mE = 1L, maxDiag = 1e4, stopIfDecreasing = TRUE, traits = "") {
    .Call(`_statgenQTLxT_EMFA`, y, k, size_param_x, cmHet, dmHet, tolerance, maxIter, size_param_cmStart, size_param_dmStart, mG, mE, maxDiag, stopIfDecreasing, traits)
}

#' Compute tYPY as in Zhou and Stephens eqn. 50.
#'
#' Compute \eqn{t(y) * P * y}, part of the log-likelihood functions from
#' equation 26 and 27 in Zhou and Stephens using equation 50. Equation 56, 57
#' and 58 are used to do the actual computations.
#'
#' It is assumed that X and Y have already been rotated by Uk, where Uk is such
#' that the kinship matrix K equals \eqn{K = Uk * Dk * t(Uk)}.\cr
#' The original X and Y are right multiplied by Uk, e.g. \code{Y <- Y * Uk}.
#' See Zhou and Stephens (2014), supplement.\cr
#' It is these rotated versions that are the input of this function.
#'
#' @inheritParams estEffsCPP
#'
#' @param size_param_x An optional c x n covariate matrix, c being the number
#' of covariates and n being the number of genotypes. c has to be at least one
#' (typically an intercept). No missing values are allowed.
#' @param vInv A n x p x p cube containing for each genotype l the
#' p x p matrix \eqn{v_l ^ {-1}} (in the notation of Zhou and Stephens).
#'
#' @return A numerical value for the \eqn{t(y) * P * y} part of the
#' log-likelihood function.
#'
#' @references Zhou, X. and Stephens, M. (2014). Efficient multivariate linear
#' mixed model algorithms for genome-wide association studies. Nature Methods,
#' February 2014, Vol. 11, p. 407–409
#'
#' @keywords internal
LLQuadFormDiagCPP <- function(y, vInv, size_param_x = NULL) {
    .Call(`_statgenQTLxT_LLQuadFormDiagCPP`, y, vInv, size_param_x)
}

#' Estimates for covariates
#'
#' Compute the estimates and standard errors for the covariates in the input
#' matrix W.
#'
#' @param y0 An n x p matrix of observed phenotypes, on p traits or environments
#' for n genotypes. No missing values are allowed.
#' @param w0 An n x c covariate matrix, c being the number of covariates and n
#' being the number of genotypes. c has to be at least one (typically an
#' intercept). No missing values are allowed.
#' @param x0 An n x ns matrix of marker scores. Neither missing values nor
#' non-segregating markers are allowed.
#' @param vg A p x p matrix of genetic covariances.
#' @param ve A p x p matrix of environmental covariances.
#' @param k An n x n genetic relatedness matrix.
#' @param returnSe Should standard errors and p-values be returned?
#' @param estCom Should the common SNP-effect model be fitted?
#' @param nCores An integer indicating the number of cores used for parallel
#' computation.
#'
#' @return A list containing the estimates, optionally the standard errors of
#' the estimates and corresponding p-values. If \code{estCom = TRUE} also
#' common SNP-effects, their standard errors and corresponding p-values and
#' the p-values for QtlxE are output.
#'
#' @references Zhou, X. and Stephens, M. (2014). Efficient multivariate linear
#' mixed model algorithms for genome-wide association studies. Nature Methods,
#' February 2014, Vol. 11, p. 407–409
#' @keywords internal
estEffsCPP <- function(y0, w0, x0, vg, ve, k, returnSe = TRUE, estCom = FALSE, nCores = NULL) {
    .Call(`_statgenQTLxT_estEffsCPP`, y0, w0, x0, vg, ve, k, returnSe, estCom, nCores)
}

getThr <- function(nCores = NULL) {
    .Call(`_statgenQTLxT_getThr`, nCores)
}

nearestPDQTLxT <- function(x, corr = FALSE, keepDiag = FALSE, do2eigen = TRUE, doSym = FALSE, doDykstra = TRUE, eigTol = 1e-6, convTol = 1e-7, posdTol = 1e-8, maxIter = 100L) {
    .Call(`_statgenQTLxT_nearestPDQTLxT`, x, corr, keepDiag, do2eigen, doSym, doDykstra, eigTol, convTol, posdTol, maxIter)
}