epiGWAS (version 1.0.2)

cond_prob: Computes the propensity scores

Description

In this function, and for each sample, we compute both propensity scores \(P(A=1| X)\) and \(P(A=0| X)\). The application of the forward algorithm on the passed hmm allows us to estimate the joint probability of (A, X), for all values of the target variant A = 0, 1, 2. The Bayes formula yields the corresponding conditional probabilities. Depending on the binarization rule, we combine them to obtain the propensity scores.

Usage

cond_prob(X, target_name, hmm, binary = FALSE, ncores = 1)

Arguments

X

genotype matrix. Make sure to assign colnames(X) beforehand.

target_name

target variant name

hmm

fitted parameters of the fastPHASE hidden Markov model. The HMM model is to be fitted with the fast_HMM function.

binary

if TRUE, the target SNP values 0 and (1,2) are respectively mapped to 0 and 1. That describes a dominant mechanism. Otherwise, if FALSE, we encode a recessive mechanism where the values 0 and 1 respectively map to (0,1) and 2.

ncores

number of threads (default 1)

Value

Two-column propensity score matrix. The first column lists the propensity score \(P\left(A=0| X\right)\), while the second gives \(P\left(A=1| X\right)\).

See Also

fast_HMM

Examples

# NOT RUN {
p <- 3 # Number of states
K <- 2 # Dimensionality of the latent space

p_init <- rep(1 / K, K)
p_trans <- array(runif((p - 1) * K * K), c(p - 1, K, K))
# Normalizing the transition probabilities
for (j in seq_len(p - 1)) {
  p_trans[j, , ] <- p_trans[j, , ] / (matrix(rowSums(p_trans[j, , ]), ncol = 1) %*% rep(1, K))
}

p_emit <- array(stats::runif(p * 3 * K), c(p, 3, K))
# Normalizing the emission probabilities
for (j in seq_len(p)) {
  p_emit[j, , ] <- p_emit[j, , ] / (matrix(rep(1, 3), ncol = 1) %*% colSums(p_emit[j, , ]))
}

hmm <- list(pInit = p_init, Q = p_trans, pEmit = p_emit)

n <- 2
X <- matrix((runif(n * p, min = 0, max = 1) < 0.4) +
            (runif(n * p, min = 0, max = 1) < 0.4),
            nrow = 2, dimnames = list(NULL, paste0("SNP_", seq_len(p))))

cond_prob(X, "SNP_2", hmm, ncores = 1, binary = TRUE)

# }