R/brm_pval.R

In Ajfrick/ajfhelpR: General Assistance for EDA and table generation

#' Generate empirical P-Value from MCMC generated by brm function
#'
#' @param brmfit Model fit by brm function
#' @return tibble with 2 columns containing Fixed Effects names and pvalues
#' @examples
#' require(brms)
#' x = 1:100
#' y = x^2 + rnorm(length(x))
#' g = rep(1:4,each = 25)
#' df = tibble(x = x,y = y,g = g)
#'
#' brmfit = brm(y ~ x^2 + (1|g), data = df)
#' brm_pval(brmfit)


brm_pval = function(fit){

  npars = length(grep("^b_", names(fit$fit@sim$samples[[1]])))-1
  samples =fit$fit@sim$samples

  if(npars>1){
    draws = rbind(as.data.frame(samples[[1]])[,2:(1+npars)],
                  as.data.frame(samples[[2]])[,2:(1+npars)],
                  as.data.frame(samples[[3]])[,2:(1+npars)],
                  as.data.frame(samples[[4]])[,2:(1+npars)])

    names = colnames(draws)
    out = tibble::tibble(names = names)
    chains = length(fit$fit@sim$samples)

    pvals = numeric(npars)

    for(i in 1:(npars)){
      coef = out$names[i]
      vals = draws %>%
        dplyr::select(which(colnames(draws)==coef))
      test = sum(vals[,1] < 0)/nrow(draws)
      pvals[i] = ifelse(median(vals[,1])>0,test,1-test)*2
    }
    out = out %>%
      dplyr::mutate(pvals = pvals)
    return(out)

  }
  if(npars == 1){

    draws = c(as.data.frame(samples[[1]])[,2:(1+npars)],
              as.data.frame(samples[[2]])[,2:(1+npars)],
              as.data.frame(samples[[3]])[,2:(1+npars)],
              as.data.frame(samples[[4]])[,2:(1+npars)])

    coef = names(samples[[1]][2])
    test = sum(draws < 0)/length(draws)
    pvals = ifelse(mean(draws)>0,test,1-test)*2
    tibble(names = coef,
           pvals = pvals)
  }
}