R/v6_mat_and_prior_funcs.R

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

v6_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
  }

  ll_mat <- ll_v6_mat(rho, b,
                      grid[,1],
                      data$beta_hat_1,
                      data$beta_hat_2,
                      data$seb1,
                      data$seb2)
  cols <- c(0, rep(1:3, each=nrow(grid)))
  if(type=="null"){
    ll_mat <- ll_mat[, cols!=2]
  }else if(type=="causal"){
    ll_mat <- ll_mat[, cols!=1]
  }
  return(ll_mat)
}

v6_prior_func <- function(pars, type=c("null", "causal", "shared"),
                          z_prior_func=NULL,
                          b_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){
      log((b/0.4)^2) + dnorm(b, 0, 0.4, 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)
  }
  return(ll)
}

v6_q_gamma_from_pi <- function(pi, type=c("null", "causal", "shared")){
  n <- length(pi)
  if(type=="null"){
    L <- (n-1)/2
    stopifnot(round(L)==L)
    gamma <- c(pi[1], sum(pi[2:(L+1)]),
               0,sum(pi[(L + 2):n]))
    q <- gamma[3]/(sum(gamma[2:3]))
    pi1 <- pi[2:(L+1)]/gamma[2]
    pi2 <- rep(NA, L)
    pi3 <- pi[(L + 2):n]/gamma[4]
    return(list("gamma"=gamma, "q"=q,
                "pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
  }
  if(type=="causal"){
    L <- (n-1)/2
    stopifnot(round(L)==L)
    gamma <- c(pi[1], 0, sum(pi[2:(L+1)]),
               sum(pi[(L + 2):n]))
    q <- gamma[3]/(sum(gamma[2:3]))
    pi1 <- rep(NA, L)
    pi2 <- pi[2:(L+1)]/gamma[3]
    pi3 <- pi[(L + 2):n]/gamma[4]
    return(list("gamma"=gamma, "q"=q,
                "pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
  }

  if(type=="shared"){
    L <- (n-1)/3
    stopifnot(round(L)==L)
    gamma <- c(pi[1], sum(pi[2:(L+1)]),
               sum(pi[(L+2):(2*L + 1)]),
               sum(pi[(2*L + 2):n]))
    q <- gamma[3]/(sum(gamma[2:3]))
    pi1 <- pi[2:(L+1)]/gamma[2]
    pi2 <- pi[(L+2):(2*L + 1)]/gamma[3]
    pi3 <- pi[(2*L + 2):n]/gamma[4]
    return(list("gamma"=gamma, "q"=q,
                "pi1"=pi1, "pi2"=pi2, "pi3"=pi3))
  }
}


v6_li_func_share <- function(pars, data, grid, gamma_red,
                          z_prior_func, b_prior_func){
  #pars: z, b, logitq
  z <- pars[1]
  rho <- tanh(z)
  b <- pars[2]
  logitq <- pars[3]
  q <- expit(pars[3])

  #grid should have first column as non-zero sigma and next three columns as pi1, pi2, pi3
  gamma <- c(gamma_red[1], (1-q)*gamma_red[2], q*gamma_red[2], gamma_red[3])
  pi <- c(gamma[1], gamma[2]*grid[,2], gamma[3]*grid[,3], gamma[4]*grid[,4])
  K <- length(pi)
  ll1 <- ll_v6(rho, b, gamma,
               grid[,2], grid[,3], grid[,4], grid[,1],
               data$beta_hat_1, data$beta_hat_2,
               data$seb1, data$seb2, data$wts)
  prior <- z_prior_func(z) + b_prior_func(b) + ddirichlet1(pi, rep(c(10, 1), c(1, K-1)))
  return(ll1 + prior)
}