R/RcppExports.R

Defines functions slcor mycor simplex_to_real real_to_simplex expit logit plog_sum_exp_mat plog_sum_exp log_sum_exp_2 log_sum_exp_mat log_sum_exp dlprior_par_dprob lprior_par dlprior_dprob lprior secalc grad_rho_m grad_deltaprime_m grad_delta_m fill_pv fill_pm obj_pbnorm_genolike prior_sigma prior_mu llike_pbnorm_genolike pbnorm_dist dmvnorm llike_geno proballgeno probgeno dmulti_double dllike_genolike_dpar dproballgenolike_dprob dprobgenolike_dprob get_prob_array get_dprobgeno_dprob_array llike_genolike proballgenolike proballgenolike_old probgenolike ddprime_dqlm dr2_dqlm dD_dqlm hessian_jointgeno llike_jointgeno em_jointgeno genolike_em simplex_proj get_Amat dDprime_dprob dr2_dprob dD_dprob dllike_geno_dpar dsimplex_to_real_dx dreal_to_simplex_dy dproballgeno_dprob dprobgeno_dprob dmulti_dprob

Documented in get_prob_array pbnorm_dist slcor

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

#' Derivative of multinomial pdf
#'
#' Gradient of multinomial pdf (NOT log-pdf) with respect to prob,
#' optionally returned on the log-scale.
#'
#' @inheritParams dmulti_double
#'
#' @author David Gerard
#'
#' @noRd
dmulti_dprob <- function(x, prob, log_p = TRUE) {
    .Call(`_ldsep_dmulti_dprob`, x, prob, log_p)
}

#' Derivative of \code{\link{probgeno}(log = TRUE)} with respect to
#' \code{prob}.
#'
#' @inheritParams probgeno
#' @param log_d A logical. Should we return the log of the derivative or not?
#'
#' @author David Gerard
#'
#' @return A vector of length 4 containing the gradient of
#'     \code{\link{probgeno}()} with respect to \code{prob}.
#'
#' @noRd
dprobgeno_dprob <- function(gA, gB, K, prob) {
    .Call(`_ldsep_dprobgeno_dprob`, gA, gB, K, prob)
}

#' Derivative of \code{\link{proballgeno}()} with respect to \code{prob}.
#'
#' @inheritParams proballgeno
#'
#' @author David Gerard
#'
#' @return A vector of length 4 containing the gradient of
#'     \code{\link{proballgeno}()} with respect to \code{prob}.
#'
#' @noRd
dproballgeno_dprob <- function(gA, gB, K, prob) {
    .Call(`_ldsep_dproballgeno_dprob`, gA, gB, K, prob)
}

#' Derivative of \code{\link{real_to_simplex}()} with respect to \code{y}.
#'
#' @param y A numeric matrix.
#'
#' @author David Gerard
#'
#' @noRd
dreal_to_simplex_dy <- function(y) {
    .Call(`_ldsep_dreal_to_simplex_dy`, y)
}

#' Derivative of \code{\link{simplex_to_real}()} with respect to \code{x}.
#'
#' @param x A simplex elements.
#'
#' @author David Gerard
#'
#' @noRd
dsimplex_to_real_dx <- function(x) {
    .Call(`_ldsep_dsimplex_to_real_dx`, x)
}

#' Derivative of \code{\link{llike_geno}()} with respect to par.
#'
#'
#' @inheritParams llike_geno
#'
#' @author David Gerard
#'
#' @noRd
dllike_geno_dpar <- function(par, gA, gB, K, alpha) {
    .Call(`_ldsep_dllike_geno_dpar`, par, gA, gB, K, alpha)
}

#' Derivative of prob[[4]] - (prob[[2]] + prob[[4]]) * (prob[[3]] + prob[[4]])
#' with respect to prob.
#'
#' These derivatives are with repsect to prob[1], prob[2], and prob[3].
#' prob[4] is defined in terms of those other three values as
#' prob[4] = 1 - sum(prob[1:3]).
#'
#' @param prob Probability vector in order of (ab, Ab, aB, AB)
#'
#' @author David Gerard
#'
#' @noRd
dD_dprob <- function(prob) {
    .Call(`_ldsep_dD_dprob`, prob)
}

#' Derivative of squared correlation with respect to prob.
#'
#' These derivatives are with repsect to prob[1], prob[2], and prob[3].
#' prob[4] is defined in terms of those other three values as
#' prob[4] = 1 - sum(prob[1:3]).
#'
#' @param prob Probability vector in order of (ab, Ab, aB, AB)
#'
#' @author David Gerard
#'
#' @noRd
dr2_dprob <- function(prob) {
    .Call(`_ldsep_dr2_dprob`, prob)
}

#' Derivative of D'with respect to prob.
#'
#' These derivatives are with repsect to prob[1], prob[2], and prob[3].
#' prob[4] is defined in terms of those other three values as
#' prob[4] = 1 - sum(prob[1:3]).
#'
#' @param prob Probability vector in order of (ab, Ab, aB, AB)
#'
#' @author David Gerard
#'
#' @noRd
dDprime_dprob <- function(prob) {
    .Call(`_ldsep_dDprime_dprob`, prob)
}

#' Fixed point iteration for \code{\link{genolike_em}()}.
#'
#'
#' @author David Gerard
#'
#' @noRd
NULL

#' Get a matrix of all possible haplotype numbers
#'
#' The numbers are returned in lexicographical order.
#'
#' @param K The ploidy of the species.
#'
#' @author David Gerard
#'
#' @noRd
get_Amat <- function(K) {
    .Call(`_ldsep_get_Amat`, K)
}

simplex_proj <- function(y) {
    .Call(`_ldsep_simplex_proj`, y)
}

#' EM algorithm to estimate haplotype frequencies
#'
#' This runs an EM algorithm to obtain the maximum likelihood estimates
#' of the haplotype frequencies for two loci when one has access
#' to genotype likelihoods.
#'
#' EM Squaring is performed via the algorithm of Ravi and Roland (2008)
#' with each iteration's projection onto the simplex using the algorithm
#' of Chen and Ye (2011). Though squaring is turned off be default because
#' it usually doesn't help.
#'
#' @param p A vector of length 4. The intialization for the
#'     haplotype frequencies.
#' @param pgA The matrix of genotype log-likelihoods for locus 1.
#'     The rows index the individuals and the columns index the genotypes.
#' @param pgB The matrix of genotype log-likelihoods for locus 2.
#'     The rows index the individuals and the columns index the genotypes.
#' @param alpha The prior sample size used in the penalty.
#' @param maxit The maximum number of EM iterations.
#' @param tol The convergence tolerance.
#' @param verbose Should we output more (\code{TRUE}) or less
#'     (\code{FALSE})?
#' @param square Should we implement squared acceleratred EM (\code{TRUE})
#'     or not (\code{FALSE})?
#'
#'
#' @references
#' \itemize{
#'   \item{Varadhan, Ravi, and Christophe Roland. "Simple and globally
#'         convergent methods for accelerating the convergence of any
#'         EM algorithm." Scandinavian Journal of Statistics 35.2 (2008):
#'         335-353.}
#'   \item{Chen, Yunmei, and Xiaojing Ye. "Projection onto a simplex."
#'         arXiv preprint arXiv:1101.6081 (2011).}
#' }
#'
#' @author David Gerard
#'
#' @noRd
genolike_em <- function(p, pgA, pgB, alpha, maxit = 500L, tol = 0.0001, verbose = FALSE, square = FALSE) {
    .Call(`_ldsep_genolike_em`, p, pgA, pgB, alpha, maxit, tol, verbose, square)
}

#' EM algorithm to estimate joint genotype frequencies.
#'
#' Implements an EM algorithm to calculate the joing distribution of dosages
#' using genotype likelihoods.
#'
#' @param p A matrix of proportions. The initialization of the joint genotype
#'     frequencies. \code{p[i,j]} is the initialization of the probability
#'     of genotype \code{i-1} on locus 1 and genotype \code{j-1} on locus 2.
#' @param pgA The matrix of genotype log-likelihoods for locus 1.
#'     The rows index the individuals and the columns index the genotypes.
#' @param pgB The matrix of genotype log-likelihoods for locus 2.
#'     The rows index the individuals and the columns index the genotypes.
#' @param alpha A matrix of prior sample sizes used as the penalty.
#' @param maxit The maximum number of EM iterations.
#' @param tol The convergence tolerance.
#' @param verbose Should we output the progress of each iteration
#'     (\code{TRUE}) or not (\code{FALSE})?
#'
#' @author David Gerard
#'
#' @noRd
#'
em_jointgeno <- function(p, pgA, pgB, alpha, maxit = 500L, tol = 0.01, verbose = FALSE) {
    .Call(`_ldsep_em_jointgeno`, p, pgA, pgB, alpha, maxit, tol, verbose)
}

#' Likelihood being maximized in \code{\link{em_jointgeno}()}
#'
#' @inheritParams em_jointgeno
#'
#' @author David Gerard
#'
#' @noRd
llike_jointgeno <- function(p, pgA, pgB, alpha) {
    .Call(`_ldsep_llike_jointgeno`, p, pgA, pgB, alpha)
}

#' Hessian of \code{\link{llike_jointgeno}()}
#'
#' Derivative of log-likelihood with respect to q_{ij} and q_{km}. The
#' ordering of the matrix is the rows index with i going fastest, and the
#' columns index with k going fastest.
#'
#' @inheritParams em_jointgeno
#'
#' @author David Gerard
#'
#' @noRd
hessian_jointgeno <- function(p, pgA, pgB, alpha) {
    .Call(`_ldsep_hessian_jointgeno`, p, pgA, pgB, alpha)
}

#' Derivative of \code{\link{Dfromg}()} with respect to gmat.
#'
#' @param p Element (i, j) is the probability of genotype i at locus 1
#'     and genotype j at locus 2.
#'
#' @author David Gerard
#'
#' @noRd
#'
dD_dqlm <- function(p) {
    .Call(`_ldsep_dD_dqlm`, p)
}

#' Gradient of squared correlation with respect to the qlm's
#'
#' @param p Element (i, j) is the probability of genotype i at locus 1
#'     and genotype j at locus 2.
#' @param dgrad The output of \code{\link{dD_dqlm}()}.
#' @param D The value of D.
#'
#' @noRd
dr2_dqlm <- function(p, dgrad, D) {
    .Call(`_ldsep_dr2_dqlm`, p, dgrad, D)
}

#' Gradient of standardized composite LD coefficient with respect to the qlm's
#'
#' @param p Element (i, j) is the probability of genotype i at locus 1
#'     and genotype j at locus 2.
#' @param dgrad The output of \code{\link{dD_dqlm}()}.
#' @param D The value of D.
#' @param Dm The value of Dm
#'
#' @noRd
ddprime_dqlm <- function(p, dgrad, D, Dm) {
    .Call(`_ldsep_ddprime_dqlm`, p, dgrad, D, Dm)
}

#' Probability of genotype likelihoods given haplotype frequencies for a
#' single individual.
#'
#' The ploidy of the species is assumed to be one less the length of
#' \code{pgA} and \code{pgB} (which must be the same length).
#'
#' @param pgA The genotype log-likelihoods at locus 1. pgA[i] is the
#'     log probability of the data given the genotype at locus 1 is i.
#' @param pgB The genotype log-likelihoods at locus 2. pgA[i] is the
#'     log probability of the data given the genotype at locus 2 is i.
#' @param prob The vector of probabilities for haplotypes (ab, Ab, aB, AB).
#' @param log_p A logical. Should we return the log probability or not?
#'
#'
#' @author David Gerard
#'
#' @noRd
probgenolike <- function(pgA, pgB, prob, log_p = TRUE) {
    .Call(`_ldsep_probgenolike`, pgA, pgB, prob, log_p)
}

#' Probability of genotype likelihoods given haplotype frequencies for all
#' individuals.
#'
#' @param pgA The matrix of genotype log-likelihoods for locus 1.
#'     The rows index the individuals and the columns index the genotypes.
#' @param pgA The matrix of genotype log-likelihoods for locus 2.
#'     The rows index the individuals and the columns index the genotypes.
#' @param prob The vector of probabilities for haplotypes (ab, Ab, aB, AB).
#' @param log_p A logical. Should we return the log probability or not?
#'
#'
#' @author David Gerard
#'
#' @noRd
proballgenolike_old <- function(pgA, pgB, prob, log_p = TRUE) {
    .Call(`_ldsep_proballgenolike_old`, pgA, pgB, prob, log_p)
}

#' Probability of genotype likelihoods given haplotype frequencies for all
#' individuals.
#'
#' @param pgA The matrix of genotype log-likelihoods for locus 1.
#'     The rows index the individuals and the columns index the genotypes.
#' @param pgA The matrix of genotype log-likelihoods for locus 2.
#'     The rows index the individuals and the columns index the genotypes.
#' @param prob The vector of probabilities for haplotypes (ab, Ab, aB, AB).
#' @param log_p A logical. Should we return the log probability or not?
#'
#' @author David Gerard
#'
#' @noRd
proballgenolike <- function(pgA, pgB, prob, log_p = TRUE) {
    .Call(`_ldsep_proballgenolike`, pgA, pgB, prob, log_p)
}

#' Likelihood function when estimating LD from genotype log-likelihoods
#'
#' @param par A vector of length 3 containing real numbers that are to
#'     be transformed into the simplex of prob (ab, Ab, aB, AB).
#' @param alpha The prior sample size used in the penalty.
#' @inheritParams proballgenolike
#'
#' @author David Gerard
#'
#' @noRd
llike_genolike <- function(par, pgA, pgB, alpha) {
    .Call(`_ldsep_llike_genolike`, par, pgA, pgB, alpha)
}

#' Obtain a matrix of derivatives of p(geno) (NOT log(p(geno)))
#' with respect to prob for all genotypes.
#'
#'
#' @param K the ploidy
#' @param prob Haplotype frequencies in order (ab, Ab, aB, AB).
#'
#' @return Element (i,j,k) is the derivative of \code{\link{probgeno}()}
#'     when gA = i, gB = j with respect to prob[k]
#'
#' @author David Gerard
#'
#' @noRd
get_dprobgeno_dprob_array <- function(K, prob) {
    .Call(`_ldsep_get_dprobgeno_dprob_array`, K, prob)
}

#' Obtain the distribution of genotypes given haplotype frequencies under HWE
#'
#' This function will calculate the (log) probabilities for all genotype
#' combinations at two loci given just the haplotype frequencies. This
#' is under the assumptions of HWE.
#'
#' @param K The ploidy of the species.
#' @param prob Haplotype frequencies in the order of (ab, Ab, aB, AB).
#' @param log_p A logical. Should we return the log-probabilities (\code{TRUE})
#'     or the probabilities (\code{FALSE}). Defaults to \code{TRUE}.
#'
#' @return Element (i, j) is the (log) probability of genotype i-1 at locus 1
#'     and genotype j-1 at locus 2.
#'
#' @author David Gerard
#'
#' @examples
#' get_prob_array(K = 6, prob = c(0.1, 0.2, 0.3, 0.4), log_p = FALSE)
#'
#' @export
#'
get_prob_array <- function(K, prob, log_p = TRUE) {
    .Call(`_ldsep_get_prob_array`, K, prob, log_p)
}

#' Gradient of \code{\link{probgenolike}()} with respect to \code{prob}.
#'
#' @inheritParams probgenolike
#'
#' @author David Gerard
#'
#' @noRd
dprobgenolike_dprob <- function(pgA, pgB, prob) {
    .Call(`_ldsep_dprobgenolike_dprob`, pgA, pgB, prob)
}

#' Gradient of \code{\link{proballgenolike}(,log = TRUE)} with respect to
#' \code{prob}
#'
#' @inheritParams proballgenolike
#'
#' @author David Gerard
#'
#' @noRd
dproballgenolike_dprob <- function(pgA, pgB, prob) {
    .Call(`_ldsep_dproballgenolike_dprob`, pgA, pgB, prob)
}

#' Derivative of \code{\link{llike_genolike}()} with respect to \code{par}.
#'
#' @inheritParams llike_genolike
#'
#' @author David Gerard
#'
#' @noRd
dllike_genolike_dpar <- function(par, pgA, pgB, alpha) {
    .Call(`_ldsep_dllike_genolike_dpar`, par, pgA, pgB, alpha)
}

#' Multinomial pdf
#'
#' @param x A vector of counts
#' @param prob A vector of probabilities
#' @param log_p A logical. Should we return the log (\code{TRUE}) or not
#'     (\code{FALSE})?
#'
#' @author David Gerard
#'
#' @noRd
dmulti_double <- function(x, prob, log_p = TRUE) {
    .Call(`_ldsep_dmulti_double`, x, prob, log_p)
}

#' Probability of genotypes given haplotype frequencies for one individual
#'
#' @param gA The number of A alleles.
#' @param gB The number of B alleles.
#' @param K The ploidy of the species.
#' @param prob A numeric vector with the probabilities of haplotypes
#'     (in order) of (ab, Ab, aB, AB).
#' @param log_p A logical. log_p A logical. Should we return the log
#'      (\code{TRUE}) or not (\code{FALSE})?
#'
#' @author David Gerard
#'
#' @noRd
probgeno <- function(gA, gB, K, prob, log_p = TRUE) {
    .Call(`_ldsep_probgeno`, gA, gB, K, prob, log_p)
}

#' Probability of genotypes given haplotype frequencies for all individuals
#'
#' @param gA Vector of number of A alleles.
#' @param gB Vector of number of B alleles.
#' @inheritParams probgeno
#'
#' @author David Gerard
#'
#' @noRd
proballgeno <- function(gA, gB, K, prob, log_p = TRUE) {
    .Call(`_ldsep_proballgeno`, gA, gB, K, prob, log_p)
}

#' Likelihood function when estimating LD directly from genotypes
#'
#' @param par A vector of length 3 containing real numbers that are to
#'     be transformed into the simplex of prob (ab, Ab, aB, AB).
#' @param alpha The prior sample size used in the penalty.
#' @inheritParams proballgeno
#'
#' @author David Gerard
#'
#' @noRd
llike_geno <- function(par, gA, gB, K, alpha) {
    .Call(`_ldsep_llike_geno`, par, gA, gB, K, alpha)
}

#' Derivative of \code{\link{dmvnorm}()} (not log) with respect to \code{R}
#' where \code{sigma = solve(R) \%*\% solve(t(R))}
#'
#'
#' @inheritParams dmvnorm
#' @param R The inverse of the cholesky square root of sigma.
#'     \code{sigma = solve(R) \%*\% solve(t(R))}
#'
#' @noRd
NULL

#' Density function of the multivariate normal distribution
#'
#' @param x A vector containing the quantile.
#' @param mu A vector containing the mean.
#' @param sigma A positive definite covariance matrix
#' @param log A logical. If \code{TRUE}, log probabilities are returned.
#'
#' @author David Gerard
#'
#' @noRd
dmvnorm <- function(x, mu, sigma, log = FALSE) {
    .Call(`_ldsep_dmvnorm`, x, mu, sigma, log)
}

#' Returns distribution of proportional bivariate normal.
#'
#' @param mu A vector of length 2 containing the mean.
#' @param sigma A 2-by-2 positive definite covariance matrix
#' @param K The ploidy of the individual.
#' @param log A logical. If \code{TRUE}, log probabilities are returned.
#'
#' @return A matrix. Element (i,j) is the (log) probability of genotype
#'     i-1 at locus 1 and j-1 at locus 2.
#'
#' @author David Gerard
#'
#' @export
#'
pbnorm_dist <- function(mu, sigma, K, log = FALSE) {
    .Call(`_ldsep_pbnorm_dist`, mu, sigma, K, log)
}

#' Log-likelihood for joint genotype distribution when using a proportional
#' normal model.
#'
#' @inheritParams pbnorm_dist
#' @param pgA The matrix of genotype log-likelihoods for locus 1.
#'     The rows index the individuals and the columns index the genotypes.
#' @param pgB The matrix of genotype log-likelihoods for locus 2.
#'     The rows index the individuals and the columns index the genotypes.
#'
#' @author David Gerard
#'
#' @noRd
llike_pbnorm_genolike <- function(pgA, pgB, mu, sigma) {
    .Call(`_ldsep_llike_pbnorm_genolike`, pgA, pgB, mu, sigma)
}

prior_mu <- function(mu, K) {
    .Call(`_ldsep_prior_mu`, mu, K)
}

prior_sigma <- function(lvec) {
    .Call(`_ldsep_prior_sigma`, lvec)
}

#' Wrapper for \code{\link{llike_pbnorm_genolike}()} that tkaes
#' a vector of parameters as input for optim().
#'
#' Also includes prior probs
#'
#' @inheritParams llike_pbnorm_genolike
#' @param par A vector of length 5. The first two elements are \code{mu}. The
#'     last three elements are c(l11, l12, l22), the lower three elements of
#'     the cholesky decomposition of sigma.
#'
#' @author David Gerard
#'
#' @noRd
#'
obj_pbnorm_genolike <- function(par, pgA, pgB) {
    .Call(`_ldsep_obj_pbnorm_genolike`, par, pgA, pgB)
}

fill_pm <- function(pm, gp) {
    invisible(.Call(`_ldsep_fill_pm`, pm, gp))
}

fill_pv <- function(pv, pm, gp) {
    invisible(.Call(`_ldsep_fill_pv`, pv, pm, gp))
}

grad_delta_m <- function(M, grad, pd) {
    invisible(.Call(`_ldsep_grad_delta_m`, M, grad, pd))
}

grad_deltaprime_m <- function(M, grad, pd) {
    invisible(.Call(`_ldsep_grad_deltaprime_m`, M, grad, pd))
}

grad_rho_m <- function(M, grad) {
    invisible(.Call(`_ldsep_grad_rho_m`, M, grad))
}

#' Calculate just the standard errors from genotype posterior array.
#'
#' Only pairwise complete observations are used to calculate standard errors.
#'
#' @param gp A three-way array with dimensions SNPs by individuals by dosage.
#'     That is, \code{gp[i, j, k]} is the posterior probability of
#'     dosage \code{k-1} for individual \code{j} at SNP \code{i}.
#' @param pm_mat The matrix of posterior mean genotypes for each individual.
#'     Rows index SNPs and columns index individuals.
#' @param pv_mat The matrix of posterior variances for each individual.
#'     Rows index SNPs and columns index individuals.
#' @param type a = D, b = r, c = D'
#'
#' @author David Gerard
#'
#' @noRd
secalc <- function(gp, pm_mat, pv_mat, type) {
    .Call(`_ldsep_secalc`, gp, pm_mat, pv_mat, type)
}

#' Prior probability for haplotype frequencies.
#'
#' @param prob The vector of probabilities for haplotypes (ab, Ab, aB, AB).
#' @param alpha The prior counts for \code{prob}.
#'
#' @author David Gerard
#'
#' @noRd
lprior <- function(prob, alpha) {
    .Call(`_ldsep_lprior`, prob, alpha)
}

#' Derivative of lprior with respect to prob.
#'
#' @param prob The vector of probabilities for haplotypes (ab, Ab, aB, AB).
#' @param alpha The prior counts for \code{prob}.
#'
#' @author David Gerard
#'
#' @noRd
dlprior_dprob <- function(prob, alpha) {
    .Call(`_ldsep_dlprior_dprob`, prob, alpha)
}

#' Prior probability on real-line scale.
#'
#' @param par A vector of length 3 containing real numbers that are to
#'     be transformed into the simplex of prob (ab, Ab, aB, AB).
#' @param alpha The prior counts for \code{prob}.
#'
#' @author David Gerard
#'
#' @noRd
lprior_par <- function(par, alpha) {
    .Call(`_ldsep_lprior_par`, par, alpha)
}

#' Derivative of \code{\link{lprior_par}()} with respect to \code{par}.
#'
#' @inheritParams lprior_par
#'
#' @author David Gerard
#'
#' @noRd
dlprior_par_dprob <- function(par, alpha) {
    .Call(`_ldsep_dlprior_par_dprob`, par, alpha)
}

#' Log-sum-exponential trick.
#'
#' @param x A vector to log-sum-exp.
#'
#' @return The log of the sum of the exponential
#'     of the elements in \code{x}.
#'
#' @author David Gerard
#'
#' @noRd
log_sum_exp <- function(x) {
    .Call(`_ldsep_log_sum_exp`, x)
}

#' log-sum-expontential of a matrix.
#'
#' @noRd
log_sum_exp_mat <- function(x) {
    .Call(`_ldsep_log_sum_exp_mat`, x)
}

#' Log-sum-exponential trick using just two doubles.
#'
#' @param x A double.
#' @param y Another double.
#'
#' @return The log of the sum of the exponential of x and y.
#'
#' @author David Gerard
#'
#' @noRd
log_sum_exp_2 <- function(x, y) {
    .Call(`_ldsep_log_sum_exp_2`, x, y)
}

#' Pair-wise log-sum-exponential
#'
#' Does pair-wise log-sum-exponential on two vectors.
#'
#' @param x A numeric vector.
#' @param y Another numeric vector.
#'
#' @author David Gerard
#'
#' @noRd
plog_sum_exp <- function(x, y) {
    .Call(`_ldsep_plog_sum_exp`, x, y)
}

#' Parallel log-sum-exp of two matrices
#'
#' @noRd
plog_sum_exp_mat <- function(x, y) {
    .Call(`_ldsep_plog_sum_exp_mat`, x, y)
}

#' The logit function.
#'
#' @param x A double between 0 and 1.
#'
#' @return The logit of \code{x}.
#'
#' @author David Gerard
#'
#' @keywords internal
#' @noRd
#'
logit <- function(x) {
    .Call(`_ldsep_logit`, x)
}

#' The expit (logistic) function.
#'
#' @param x A double.
#'
#' @return The expit (logistic) of \code{x}.
#'
#' @keywords internal
#' @noRd
#'
#' @author David Gerard
expit <- function(x) {
    .Call(`_ldsep_expit`, x)
}

#' Convert real line to simplex using Stan technique
#'
#' @param y A vector of numbers of length K-1
#'
#' @return A vector on the simplex of length K
#'
#' @author David Gerard
#'
#' @seealso \code{\link{simplex_to_real}()} for the inverse function.
#'
#' @noRd
real_to_simplex <- function(y) {
    .Call(`_ldsep_real_to_simplex`, y)
}

#' Convert simplex to real-line using Stan technique
#'
#' @param x A vector on the simplex of length K.
#'
#' @return A vector of numbers of length K-1
#'
#' @author David Gerard
#'
#' @seealso \code{\link{real_to_simplex}()} for the inverse function.
#'
#' @noRd
simplex_to_real <- function(x) {
    .Call(`_ldsep_simplex_to_real`, x)
}

#' Pearson correlation between x and y using pairwise complete obs
#'
#' @noRd
mycor <- function(x, y) {
    .Call(`_ldsep_mycor`, x, y)
}

#' Sliding window correlation
#'
#' Calculates the pairwise Pearson correlation between all columns
#' within a fixed window size (\code{win})
#' using the \code{use = "pairwise.complete.obs"} option
#' from \code{\link[stats]{cor}()}. That is, the correlation
#' between each pair of variables is computed using all complete pairs
#' of observations on those variables.
#'
#' @param x A numeric matrix. The variables index the columns.
#' @param win The size of the window. Defaults to 1.
#'
#' @return A correlation matrix with only the observations within a window
#'     containing calculated correlations.
#'
#' @export
#'
#' @examples
#' set.seed(1)
#' n <- 10
#' p <- 100
#' xmat <- matrix(rnorm(n * p), ncol = n)
#' xmat[sample(n * p, size = 30)] <- NA_real_
#' slcor(xmat, win = 2)
#'
#' @author David Gerard
slcor <- function(x, win = 1L) {
    .Call(`_ldsep_slcor`, x, win)
}

Try the ldsep package in your browser

Any scripts or data that you put into this service are public.

ldsep documentation built on Oct. 19, 2022, 1:08 a.m.