# R/hbmr.R In BMRV: Bayesian Models for Rare Variant Association Analysis

#### Documented in hbmr

```hbmr <-
function(pheno, geno, qi= matrix(), fam = 0, kin = matrix(), iter = 10000, burnin = 500, gq = 20, imp = 0.1, cov = matrix(), maf = c(), rvinfo=FALSE, pa = 1.3, pb = 0.04)
{

if ( ( ! is.matrix(geno) ) | ( ! is.matrix(qi) ) | ( ! is.matrix(kin) ) | ( ! is.matrix(cov) )) {
stop("Genotype, quality, relatedness and covariate data must be matrices!")
}

num_sub <- length(pheno)
num_rv <- ncol(geno)
num_cov <- 0

burnin <- ceiling(burnin)
iter <- ceiling(iter)

hbmr_check(geno, qi, fam, kin, iter, burnin, imp, num_sub, num_rv, num_cov)

if(sum(dim(cov))>2)
{
if(nrow(cov)!=num_sub)
{
stop("Covariate matrix has incorrect rows!")
}
num_cov <- ncol(cov)
}
else
{
cov <- matrix(0)
}

G <- matrix(0)
if(fam == 1)
{
if(sum(dim(kin)==c(num_sub,num_sub))<2)
{
stop("The dimension of the kinship matrix is incorrect!")
}
ei_ran <- eigen(kin)
G <- ei_ran\$vectors%*%diag(sqrt(ifelse(ei_ran\$value<0,0,ei_ran\$value)))
}

if((sum(dim(qi))==2) & (is.na(qi[1,1])))
{
qi <- matrix(99, num_sub, num_rv)
}

if(length(maf)==0)
{
maf <- -1
}
else
{
if(length(maf)!=num_rv)
{
stop("The MAF information is incomplete!")
}
}

y <- pheno
x <- t(geno)
q <- t(qi)
c <- t(cov)
# ran <- matrix(0)
ran <- t(G)

rv_detail <- 0
if(rvinfo == TRUE)
{rv_detail <- 1}
if(rv_detail == 0)
{
multResult <- rep(0,6+num_cov)
}
else
{
multResult <- rep(0,6+num_cov+2*num_rv)
}

output =.C("CWrapper_hbmr",
product = as.double(multResult),
nRows = as.integer(num_sub),
nCols = as.integer(num_rv),
nCols2 = as.integer(num_cov),
fam = as.integer(fam),
matrix1 = as.double(y),
matrix2 = as.double(x),
matrix3 = as.double(q),
matrix4 = as.double(c),
kin_m = as.double(ran),
arg_m = as.double(maf),
arg_i = as.double(imp),
arg_n = as.integer(iter),
arg_bu = as.integer(burnin),
arg_t = as.double(gq),
gamma_est = as.integer(rv_detail),
arg_a = as.double(pa),
arg_b = as.double(pb)
)

#return(output\$product)

if(output\$product[1]!=1)
{
bf <- output\$product[1]/(1-output\$product[1])
}
else
{
bf <- Inf
}

if(output\$product[2]!=1)
{
bf_rb <- output\$product[2]/(1-output\$product[2])
}
else
{
bf_rb <- iter-burnin
warning("The estimated Bayes factor bf_rb is the lower bound due to the limitation of numerical precision.")
}

re_cov <- c()
if(num_cov>0)
{re_cov <- output\$product[7:(6+num_cov)]}

re_rv_es <- c()
if(rvinfo==TRUE)
{
re_rv_es <- output\$product[(6+num_cov+1):(6+num_cov+num_rv)]
}

re_rv_sd <- c()
if(rvinfo==TRUE)
{
re_rv_sd <- output\$product[(6+num_cov+num_rv+1):(6+num_cov+num_rv+num_rv)]
}

p_val <- NA
if(bf_rb>2)
{
fun <- function (x) (-1)/(exp(1)*x*log(x))-bf_rb
p_val <- uniroot(fun, tol = 0.00000000000000000001, interval=c(0.0000000000000000001,0.5))
}

return(list(BF = bf, BF_RB = bf_rb, p_upper = p_val['root'], mean = output\$product[3], var=output\$product[4], est_geno = output\$product[5],
var_ran = output\$product[6], rv_mean_es = re_rv_es, rv_sd_es = re_rv_sd, mean_cov=re_cov))

}
```

## Try the BMRV package in your browser

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

BMRV documentation built on May 2, 2019, 3:28 a.m.