R/RcppExports.R

In pedmod: Pedigree Models

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

.biconnected_components <- function(from, to, weights_ids, weights, edge_weights) {
    .Call(`_pedmod_get_biconnected_components`, from, to, weights_ids, weights, edge_weights)
}

.block_cut_tree <- function(from, to, weights_ids, weights, edge_weights) {
    .Call(`_pedmod_get_block_cut_tree`, from, to, weights_ids, weights, edge_weights)
}

.max_balanced_partition <- function(from, to, weights_ids, weights, edge_weights, slack, max_kl_it_inner, max_kl_it, trace, check_weights, do_reorder) {
    .Call(`_pedmod_get_max_balanced_partition`, from, to, weights_ids, weights, edge_weights, slack, max_kl_it_inner, max_kl_it, trace, check_weights, do_reorder)
}

.unconnected_partition <- function(from, to, weights_ids, weights, edge_weights, slack, max_kl_it_inner, max_kl_it, trace, init) {
    .Call(`_pedmod_unconnected_partition_rcpp`, from, to, weights_ids, weights, edge_weights, slack, max_kl_it_inner, max_kl_it, trace, init)
}

#' Multivariate Normal Distribution CDF and Its Derivative
#' @description
#' Provides an approximation of the multivariate normal distribution CDF
#' over a hyperrectangle and the derivative with respect to the mean vector
#' and the covariance matrix.
#'
#' @param lower numeric vector with lower bounds.
#' @param upper numeric vector with upper bounds.
#' @param mu numeric vector with means.
#' @param sigma covariance matrix.
#' @param maxvls maximum number of samples in the approximation.
#' @param abs_eps absolute convergence threshold.
#' @param rel_eps relative convergence threshold.
#' @param minvls minimum number of samples. Negative values provides a
#' default which depends on the dimension of the integration.
#' @param do_reorder \code{TRUE} if a heuristic variable reordering should
#' be used. \code{TRUE} is likely the best value.
#' @param use_aprx \code{TRUE} if a less precise approximation of
#' \code{\link{pnorm}} and \code{\link{qnorm}} should be used. This may
#' reduce the computation time while not affecting the result much.
#' @param method integer with the method to use. Zero yields randomized Korobov
#' lattice rules while one yields scrambled Sobol sequences.
#' @param n_sequences number of randomized quasi-Monte Carlo sequences to use.
#' More samples yields a better estimate of the error but a worse
#' approximation. Eight is used in the original Fortran code. If one is
#' used then the error will be set to zero because it cannot be estimated.
#' @param use_tilting \code{TRUE} if the minimax tilting method suggested
#' by Botev (2017) should be used. See \doi{10.1111/rssb.12162}.
#'
#' @return
#' \code{mvndst:}
#' An approximation of the CDF. The \code{"n_it"} attribute shows the number of
#' integrand evaluations, the \code{"inform"} attribute is zero if the
#' requested precision is achieved, and the \code{"abserr"} attribute
#' shows 3.5 times the estimated standard error.
#'
#' @examples
#' # simulate covariance matrix and the upper bound
#' set.seed(1)
#' n <- 10L
#' S <- drop(rWishart(1L, 2 * n, diag(n) / 2 / n))
#' u <- rnorm(n)
#'
#' system.time(pedmod_res <- mvndst(
#'     lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n),
#'     maxvls = 1e6, abs_eps = 0, rel_eps = 1e-4, use_aprx = TRUE))
#' pedmod_res
#'
#' # compare with mvtnorm
#' if(require(mvtnorm)){
#'     mvtnorm_time <- system.time(mvtnorm_res <- mvtnorm::pmvnorm(
#'         upper = u, sigma = S, algorithm = GenzBretz(
#'             maxpts = 1e6, abseps = 0, releps = 1e-4)))
#'     cat("mvtnorm_res:\n")
#'     print(mvtnorm_res)
#'
#'     cat("mvtnorm_time:\n")
#'     print(mvtnorm_time)
#' }
#'
#' # with titling
#' system.time(pedmod_res <- mvndst(
#'     lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n),
#'     maxvls = 1e6, abs_eps = 0, rel_eps = 1e-4, use_tilting = TRUE))
#' pedmod_res
#'
#' # compare with TruncatedNormal
#' if(require(TruncatedNormal)){
#'     TruncatedNormal_time <- system.time(
#'         TruncatedNormal_res <- TruncatedNormal::pmvnorm(
#'             lb = rep(-Inf, n), ub = u, sigma = S,
#'             B = attr(pedmod_res, "n_it"), type = "qmc"))
#'     cat("TruncatedNormal_res:\n")
#'     print(TruncatedNormal_res)
#'
#'     cat("TruncatedNormal_time:\n")
#'     print(TruncatedNormal_time)
#' }
#'
#' # check the gradient
#' system.time(pedmod_res <- mvndst_grad(
#'   lower = rep(-Inf, n), upper = u, sigma = S, mu = numeric(n),
#'   maxvls = 1e5, minvls = 1e5, abs_eps = 0, rel_eps = 1e-4, use_aprx = TRUE))
#' pedmod_res
#'
#' \donttest{# compare with numerical differentiation. Should give the same up to Monte
#' # Carlo and finite difference error
#' if(require(numDeriv)){
#'   num_res <- grad(
#'     function(par){
#'       set.seed(1)
#'       mu <- head(par, n)
#'       S[upper.tri(S, TRUE)] <- tail(par, -n)
#'       S[lower.tri(S)] <- t(S)[lower.tri(S)]
#'       mvndst(
#'         lower = rep(-Inf, n), upper = u, sigma = S, mu = mu,
#'         maxvls = 1e4, minvls = 1e4, abs_eps = 0, rel_eps = 1e-4,
#'         use_aprx = TRUE)
#'     }, c(numeric(n), S[upper.tri(S, TRUE)]),
#'     method.args = list(d = .01, r = 2))
#'
#'   d_mu <- head(num_res, n)
#'   d_sigma <- matrix(0, n, n)
#'   d_sigma[upper.tri(d_sigma, TRUE)] <- tail(num_res, -n)
#'   d_sigma[upper.tri(d_sigma)] <- d_sigma[upper.tri(d_sigma)] / 2
#'   d_sigma[lower.tri(d_sigma)] <- t(d_sigma)[lower.tri(d_sigma)]
#'
#'   cat("numerical derivatives\n")
#'   print(rbind(numDeriv = d_mu,
#'               pedmod = pedmod_res$d_mu))
#'   print(d_sigma)
#'   cat("\nd_sigma from pedmod\n")
#'   print(pedmod_res$d_sigma) # for comparison
#' }}
#'
#' @export
mvndst <- function(lower, upper, mu, sigma, maxvls = 25000L, abs_eps = .001, rel_eps = 0L, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, method = 0L, n_sequences = 8L, use_tilting = FALSE) {
    .Call(`_pedmod_mvndst`, lower, upper, mu, sigma, maxvls, abs_eps, rel_eps, minvls, do_reorder, use_aprx, method, n_sequences, use_tilting)
}

#' @rdname mvndst
#' @return
#' \code{mvndst_grad:}
#' A list with
#' \itemize{
#'   \item \code{likelihood}: the likelihood approximation.
#'   \item \code{d_mu}: the derivative with respect to the the mean vector.
#'   \item \code{d_sigma}: the derivative with respect to the covariance matrix
#'   ignoring the symmetry (i.e. working the \eqn{n^2} parameters with
#'   \eqn{n} being the dimension rather than the \eqn{n(n + 1) / 2}
#'   free parameters).
#' }
#'
#' @export
mvndst_grad <- function(lower, upper, mu, sigma, maxvls = 25000L, abs_eps = .001, rel_eps = 0L, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, method = 0L, n_sequences = 8L, use_tilting = FALSE) {
    .Call(`_pedmod_mvndst_grad`, lower, upper, mu, sigma, maxvls, abs_eps, rel_eps, minvls, do_reorder, use_aprx, method, n_sequences, use_tilting)
}

#' Get a C++ Object for Log Marginal Likelihood Approximations
#'
#' @description
#' Constructs an object needed for \code{\link{eval_pedigree_ll}} and
#' \code{\link{eval_pedigree_grad}}.
#'
#' @param data \code{\link{list}} where each element is a list for a cluster
#' with an:
#' \itemize{
#'   \item{\code{"X"}}{ element with the design matrix for the fixed effect,}
#'   \item{\code{"Z"}}{ element with the design matrix for the loadings of the effects (only needed for \code{pedigree_ll_terms_loadings}),}
#'   \item{\code{"y"}}{ element with the zero-one outcomes, and}
#'   \item{\code{"scale_mats"}}{ element with a list where each element is a
#' scale/correlation matrix for a particular type of effect.}
#' }
#' @param max_threads maximum number of threads to use.
#' @param n_sequences number of randomized quasi-Monte Carlo sequences to use.
#' More samples yields a better estimate of the error but a worse
#' approximation. Eight is used in the original Fortran code. If one is
#' used then the error will be set to zero because it cannot be estimated.
#'
#' @details
#' An intercept column is not added to the \code{X} matrices
#' like what \code{\link{lm.fit}} and \code{\link{glm.fit}} do.
#' Thus, it is often important that the user adds an intercept column
#' to these matrices as it is hardly ever justified to not include the
#' intercept (the exceptions being e.g. when splines are used which include
#' the intercept and with certain dummy designs). This equally holds for
#' the \code{Z} matrices with \code{pedigree_ll_terms_loadings}.
#'
#' \code{pedigree_ll_terms_loadings} relax the assumption that the scale
#' parameter is the same for all individuals. \code{pedigree_ll_terms_loadings}
#' and \code{pedigree_ll_terms} yield the same model if \code{"Z"} is an
#' intercept column for all families but with a different parameterization.
#' In this case, \code{pedigree_ll_terms} will be
#' faster. See \code{vignette("pedmod", "pedmod")} for examples of using
#' \code{pedigree_ll_terms_loadings}.
#'
#' @examples
#' # three families as an example
#' fam_dat <- list(
#'   list(
#'     y = c(FALSE, TRUE, FALSE, FALSE),
#'     X = structure(c(
#'       1, 1, 1, 1, 1.2922654151273, 0.358134905909256, -0.734963997107464,
#'       0.855235473516044, -1.16189500386223, -0.387298334620742,
#'       0.387298334620742, 1.16189500386223),
#'       .Dim = 4:3, .Dimnames = list( NULL, c("(Intercept)", "X1", ""))),
#'     rel_mat = structure(c(
#'       1, 0.5, 0.5, 0.125, 0.5, 1, 0.5, 0.125, 0.5, 0.5,
#'       1, 0.125, 0.125, 0.125, 0.125, 1), .Dim = c(4L, 4L)),
#'     met_mat = structure(c(1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1),
#'                         .Dim = c(4L, 4L))),
#'   list(
#'     y = c(FALSE, FALSE, FALSE),
#'     X = structure(c(
#'       1, 1, 1, -0.0388728997202442, -0.0913782435233639,
#'       -0.0801619722392612, -1, 0, 1), .Dim = c(3L, 3L)),
#'     rel_mat = structure(c(
#'       1, 0.5, 0.125, 0.5, 1, 0.125, 0.125, 0.125, 1), .Dim = c(3L, 3L)),
#'     met_mat = structure(c(
#'       1, 1, 0, 1, 1, 0, 0, 0, 1), .Dim = c(3L, 3L))),
#'   list(
#'     y = c(TRUE, FALSE),
#'     X = structure(c(
#'       1, 1, 0.305275750370738, -1.49482995913648,  -0.707106781186547,
#'       0.707106781186547),
#'       .Dim = 2:3, .Dimnames = list( NULL, c("(Intercept)", "X1", ""))),
#'     rel_mat = structure(c(1, 0.5,  0.5, 1), .Dim = c(2L, 2L)),
#'     met_mat = structure(c(1, 1, 1, 1), .Dim = c(2L,  2L))))
#'
#' # get the data into the format needed for the package
#' dat_arg <- lapply(fam_dat, function(x){
#'   # we need the following for each family:
#'   #   y: the zero-one outcomes.
#'   #   X: the design matrix for the fixed effects.
#'   #   scale_mats: list with the scale matrices for each type of effect.
#'   list(y = as.numeric(x$y), X = x$X,
#'        scale_mats = list(x$rel_mat, x$met_mat))
#' })
#'
#' # get a pointer to the C++ object
#' ptr <- pedigree_ll_terms(dat_arg, max_threads = 1L)
#'
#' # get the argument for a the version with loadings
#' dat_arg_loadings <- lapply(fam_dat, function(x){
#'   list(y = as.numeric(x$y), X = x$X, Z = x$X[, 1:2],
#'        scale_mats = list(x$rel_mat, x$met_mat))
#' })
#'
#' ptr <- pedigree_ll_terms_loadings(dat_arg_loadings, max_threads = 1L)
#'
#' @export
pedigree_ll_terms <- function(data, max_threads = 1L, n_sequences = 8L) {
    .Call(`_pedmod_pedigree_ll_terms`, data, max_threads, n_sequences)
}

.get_n_scales <- function(ptr) {
    .Call(`_pedmod_get_n_scales`, ptr)
}

.get_n_scales_loadings <- function(ptr) {
    .Call(`_pedmod_get_n_scales_loadings`, ptr)
}

.get_n_terms <- function(ptr) {
    .Call(`_pedmod_get_n_terms`, ptr)
}

.get_n_terms_loadings <- function(ptr) {
    .Call(`_pedmod_get_n_terms_loadings`, ptr)
}

eval_pedigree_ll_cpp <- function(ptr, par, maxvls, abs_eps, rel_eps, indices = NULL, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, n_threads = 1L, cluster_weights = NULL, method = 0L, use_tilting = FALSE, vls_scales = NULL) {
    .Call(`_pedmod_eval_pedigree_ll`, ptr, par, maxvls, abs_eps, rel_eps, indices, minvls, do_reorder, use_aprx, n_threads, cluster_weights, method, use_tilting, vls_scales)
}

eval_pedigree_grad_cpp <- function(ptr, par, maxvls, abs_eps, rel_eps, indices = NULL, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, n_threads = 1L, cluster_weights = NULL, method = 0L, use_tilting = FALSE, vls_scales = NULL) {
    .Call(`_pedmod_eval_pedigree_grad`, ptr, par, maxvls, abs_eps, rel_eps, indices, minvls, do_reorder, use_aprx, n_threads, cluster_weights, method, use_tilting, vls_scales)
}

eval_pedigree_hess_cpp <- function(ptr, par, maxvls, abs_eps, rel_eps, indices = NULL, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, n_threads = 1L, cluster_weights = NULL, method = 0L, use_tilting = FALSE, vls_scales = NULL) {
    .Call(`_pedmod_eval_pedigree_hess`, ptr, par, maxvls, abs_eps, rel_eps, indices, minvls, do_reorder, use_aprx, n_threads, cluster_weights, method, use_tilting, vls_scales)
}

#' @rdname pedigree_ll_terms
#'
#' @export
pedigree_ll_terms_loadings <- function(data, max_threads = 1L, n_sequences = 8L) {
    .Call(`_pedmod_pedigree_ll_terms_loadings`, data, max_threads, n_sequences)
}

eval_pedigree_ll_loadings_cpp <- function(ptr, par, maxvls, abs_eps, rel_eps, indices = NULL, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, n_threads = 1L, cluster_weights = NULL, method = 0L, use_tilting = FALSE, vls_scales = NULL) {
    .Call(`_pedmod_eval_pedigree_ll_loadings`, ptr, par, maxvls, abs_eps, rel_eps, indices, minvls, do_reorder, use_aprx, n_threads, cluster_weights, method, use_tilting, vls_scales)
}

eval_pedigree_grad_loadings_cpp <- function(ptr, par, maxvls, abs_eps, rel_eps, indices = NULL, minvls = -1L, do_reorder = TRUE, use_aprx = FALSE, n_threads = 1L, cluster_weights = NULL, method = 0L, use_tilting = FALSE, vls_scales = NULL) {
    .Call(`_pedmod_eval_pedigree_grad_loadings`, ptr, par, maxvls, abs_eps, rel_eps, indices, minvls, do_reorder, use_aprx, n_threads, cluster_weights, method, use_tilting, vls_scales)
}