R/cor_functions.R

In abnormally-distributed/abdisttools: Random Functions for Personal Use

multivariate_means = function(data, chains = 4, cores = 4, iter = 10000, warmup = 2000, control = 
                                list(adapt_delta = .85, max_treedepth = 20)){
  
  N = nrow(data)
  P = ncol(data)
  
  infmatrix = function(data){ 
    N = nrow(data)
    info.matrix = corpcor::pseudoinverse(crossprod(as.matrix(data)))
    return(info.matrix * N) 
  }
  
  tau = infmatrix(data)
  avg.mean = mean(colMeans(data))
  prior_mean_vector = rep(avg.mean, P)
  scale = corpcor::make.positive.definite(cov(data))
  
  standata = list(N = N, P = P, y = data, prior_mean_vector = prior_mean_vector, tau = tau, scale = scale)
  
  sampling(model, standata,
           pars = c("means", "sigma", "nu", "nu_alpha","nu_beta","prior_means"),
           cores = cores, chains = chains, iter = iter, warmup = warmup, control = control)
}

t_test_bayes = function(formula, data, chains = 4, cores = 4, iter = 10000, warmup = 2000, control = 
                          list(adapt_delta = .85, max_treedepth = 20)){
  
  
  X = model.frame(formula, data)[,2]
  y = model.frame(formula, data)[,1]
  
  grp1 = which(X==1)
  grp2 = which(X==2)
  y1 = y[grp1]
  y2 = y[grp2]
  data = cbind.data.frame(Group1 = y1, Group2 = y2)
  N = nrow(data)
  P = ncol(data)
  
  infmatrix = function(data){ 
    N = nrow(data)
    info.matrix = corpcor::pseudoinverse(crossprod(as.matrix(data)))
    return(info.matrix * N) 
  }
  
  tau = infmatrix(data)
  avg.mean = mean(apply(data, 2, median))
  prior_mean_vector = rep(avg.mean, P)
  scale = corpcor::make.positive.definite(cov(data))
  
  standata = list(N = N, P = P, y = data, prior_mean_vector = prior_mean_vector, tau = tau, scale = scale)
  
  sampling(model, standata,
           pars = c("mu", "sigma", "muDiff", "CLES", "nu", "nu_alpha", "nu_beta", "prior_mu", "prior_muDiff"),
           cores = cores, chains = chains, iter = iter, warmup = warmup, control = control)
}


extract_cor_post = function(fit, ...){
  post = extract_post(fit, ...)
  wch = grep("cormat", x = colnames(post))
  N = attr(fit, "P")
  mat = matrix(0, N, N)
  cors = post[,wch]
  cors = cors[,which(lower.tri(mat) == TRUE)]
  others = post[,-wch]
  post = cbind.data.frame(cors, others)
  return(post)
}

extract_cor_prior = function(fit, ...){
  prior = extract_prior(fit, ...)
  wch = grep("prior", x = colnames(prior))
  N = attr(fit, "P")
  mat = matrix(0, N, N)
  cors = prior[,wch]
  cors = cors[,which(lower.tri(mat) == TRUE)]
  return(cors)
}

cor_summary = function(fit, digits = 2,  bayes.factor = TRUE, estimate.method = "median", cred.method = "ETI", cred.level = .95, ...) {
  post = extract_post(fit)
  prior = extract_prior(fit)
  wch = grep("cormat", x = colnames(post))
  N = attr(fit, "P")
  mat = matrix(0, N, N)
  cors = post[,wch]
  cors = cors[,which(lower.tri(mat) == TRUE)]
  prior = prior[,which(lower.tri(mat) == TRUE)]
  sum = post_summary(cors, estimate.method = estimate.method, cred.method = cred.method, cred.level = cred.level, ...)
  sum = sum %>% mutate_if(is.numeric, signif, digits)
  if (bayes.factor == TRUE) {
    BF = as.data.frame(t(sapply(1:ncol(prior), function(x) abdisttools:::hypothesis_test(cors[, x], prior[, x], H0 = 0))))
    BF = BF %>% mutate_if(is.numeric, signif, 4)
    
    sum = cbind.data.frame(sum, BF)
    sum = sum %>% transmute(term = term, estimate = estimate, std.error = std.error, cred.low = cred.low,
                            cred.high = cred.high, BF01 = BF01)
    sum = as_tibble(sum)
  }
  
  return(sum)
}

extract_cormat = function(fit) {
  P = attr(fit, "P")
  post = extract_post(fit)
  prior = extract_prior(fit)
  wch = grep("cormat", x = colnames(post))
  post = post[,wch]
  estimates = colMeans(post)
  matrix(estimates, P, P)
}