#' 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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.