R/v4_mat_and_prior_funcs.R

In jean997/sherlockAsh: CAUSE: Causal Analysis Using Summary Effect Estimates

#Returns a matrix that is number of SNPs by number of grid points
# ll_ij = p(beta_hat_1_i, beta_hat_2_i | rho, b, q, sigma_1j, sigma_2j)
v4_li_mat_func <- function(pars, data, grid, type=c("null", "causal", "shared"), ...){
  type <- match.arg(type)
  rho <- pars[1]
  if(type %in% c("causal", "shared")){
    stopifnot(length(pars) > 1)
    b <- pars[2]
  }else{
    b <- 0
  }
  if(type=="shared"){
    stopifnot(length(pars)==3)
    q <- pars[3]
  }else{
    q <- 1
  }
  ll_mat <- ll_v4_mat(rho, b, q,
                      grid$S1, grid$S2,
                      data$beta_hat_1, data$beta_hat_2,
                      data$seb1, data$seb2)
}

#Calculates p(rho, b, q)
v4_prior_func <- function(pars, type=c("null", "causal", "shared"),
                        z_prior_func=NULL, b_prior_func = NULL, q_prior_func  = NULL){

  type <- match.arg(type)
  arctanh <- function(rho){
    0.5*log((1+rho)/(1-rho))
  }

  if(is.null(z_prior_func)){
    z_prior_func = function(z){dnorm(z, 0, 0.5, log=TRUE)}
  }
  if(is.null(b_prior_func)){
    b_prior_func <- function(b){dnorm(b, 0, 0.6, log=TRUE)}
  }
  if(is.null(q_prior_func)){
    q_prior_func <- function(q){dnorm(q, 0, 1.2, log=TRUE)}
  }
  rho <- pars[1]
  z <- arctanh(rho)
  ll <- z_prior_func(z)
  if(type %in% c("causal", "shared")){
    stopifnot(length(pars) > 1)
    b <- pars[2]
    ll <- ll + b_prior_func(b)
  }else{
    b <- 0
  }
  if(type=="shared"){
    stopifnot(length(pars)==3)
    q <- pars[3]
    ll <- ll + q_prior_func(q)
  }else{
    q <- 0
  }
  return(ll)
}