Installation

remotes::install_github("crsh/jab")

Example

library("jab")library("ggplot2")# Fit regression modeldata(attitude)attitude_z<-data.frame(scale(attitude))attitude_lm<- lm(rating~0+.,data=attitude_z)attitude_tidy_lm<-broom::tidy(attitude_lm)attitude_tidy_lm#> # A tibble: 6 × 5#>   term       estimate std.error statistic  p.value#>   <chr>         <dbl>     <dbl>     <dbl>    <dbl>#> 1 complaints   0.671      0.172     3.89  0.000694#> 2 privileges  -0.0734     0.134    -0.550 0.588#> 3 learning     0.309      0.159     1.94  0.0640#> 4 raises       0.0698     0.185     0.377 0.710#> 5 critical     0.0312     0.117     0.267 0.792#> 6 advance     -0.183      0.147    -1.24  0.225# Specify prior distribution and approximate Bayes factorattitude_jab<- jab(attitude_lm  ,prior=dcauchy  ,location=0  ,scale= sqrt(2)/4)attitude_jab#> complaints privileges   learning     raises   critical    advance#> 0.03751703 2.98754241 0.88364001 2.31936704 3.68709775 1.82951234

# Calculate JZS-Bayes factorattitude_jzs<-BayesFactor::regressionBF(rating~.  ,data=attitude  ,rscaleCont= sqrt(2)/4  ,whichModels="top"  ,progress=FALSE)# Compare resultstibble::tibble(predictor=attitude_tidy_lm$term# Frequentist p-values  ,p=attitude_tidy_lm$p.value# Bayes factors in favor of the null hypothesis  ,jab=attitude_jab  ,jzs= rev(as.vector(attitude_jzs))# Naive posterior probabilities  ,jab_pp=jab/ (jab+1)  ,jzs_pp=jzs/ (jzs+1))#> # A tibble: 6 × 6#>   predictor         p    jab    jzs jab_pp jzs_pp#>   <chr>         <dbl>  <dbl>  <dbl>  <dbl>  <dbl>#> 1 complaints 0.000694 0.0375 0.0231 0.0362 0.0225#> 2 privileges 0.588    2.99   2.92   0.749  0.745#> 3 learning   0.0640   0.884  0.727  0.469  0.421#> 4 raises     0.710    2.32   3.13   0.699  0.758#> 5 critical   0.792    3.69   3.23   0.787  0.764#> 6 advance    0.225    1.83   1.73   0.647  0.634

Varying prior distributions

jab(attitude_lm  ,prior=dcauchy  ,location=0  ,scale= c(rep(0.5,3), rep(sqrt(2)/4,3)))#> complaints privileges   learning     raises   critical    advance#>  0.0322969  4.1376723  0.9791980  2.3193670  3.6870978  1.8295123

Prior sensitivity

# Specify designjab_sensitivity<- expand.grid(coef= names(coef(attitude_lm))  ,r= seq(0.2,1.5,length.out=50))|># Calculate Bayes factors for each prior settingdplyr::group_by(r)|>dplyr::mutate(jab= jab(attitude_lm      ,prior=dcauchy      ,location=0      ,scale=r    )  )# Plot resultsggplot(jab_sensitivity)+  aes(x=r,y=jab/ (1+jab),color=coef)+  geom_hline(yintercept=0.5    ,linetype="22"    ,color= grey(0.7)  )+  geom_line(linewidth=1.5)+  scale_color_viridis_d()+  lims(y= c(0,1))+  labs(x= bquote(italic(r))    ,y="Naive posterior probability"    ,color="Coefficient"  )+papaja::theme_apa(box=TRUE)

Sequential analyses

# Specify designsequential_jab<- expand.grid(coef= names(coef(attitude_lm))  ,n=10:nrow(attitude_z))|># Calculate Bayes factors for each subsampledplyr::group_by(n)|>dplyr::mutate(jab= jab(      update(attitude_lm,data=attitude_z[1:unique(n), ])      ,dcauchy      ,location=0      ,scale= sqrt(2)/4    )    ,jab_pp=jab/ (jab+1)  )# Plot resultsggplot(sequential_jab)+  aes(x=n,y=jab_pp,color=coef)+  geom_line(linewidth=1.5)+  scale_color_viridis_d()+  lims(y= c(0,1))+  labs(x= bquote(italic(n))    ,y="Naive posterior probability"    ,color="Coefficient"  )+papaja::theme_apa(box=TRUE)

What’s in a p-value?

library("geomtextpath")p_boundaries<- c(0.0001,0.001,0.01,0.05,0.1,1)dat<- expand.grid(p= exp(seq(log(0.00005), log(1),length.out=100))  ,n= exp(seq(log(3), log(10000),length.out=100)))|>  transform(jab_p=1/jab::jab_p(p,n))evidence_labels<-data.frame(n= c(17,50,75,150,350,800,2000,4800)  ,p= c(0.0002,0.0009,0.00225,0.005,0.019,0.065,0.175,0.45)  ,label= c("Extreme","Very strong","Strong","Moderate","Anecdotal","Moderate","Strong","Very strong")  ,angle=-c(17,17,17,17,17,18,21,24)+3)plot_settings<-list(  scale_x_continuous(expand= expansion(0,0)    ,breaks= c(5,10,20,50,100,250,500,1000,2500,5000,10000)    ,trans="log"    ,name= bquote("Sample size"~ italic(n))  )  , scale_y_continuous(expand= expansion(0,0)    ,breaks=p_boundaries    ,labels= format(p_boundaries,scientific=FALSE,drop0trailing=TRUE)    ,trans="log"    ,name= bquote(italic(p)*"-value")  )  , scale_fill_viridis_c(guide="none")  , theme_minimal(base_size=16)  , theme(axis.ticks.length= unit(5,"pt")    ,axis.ticks.x= element_line()    ,axis.ticks.y= element_line()    ,plot.margin= margin(0.5,0.5,0.5,0.5,"cm")    ,axis.title.x= element_text(margin= margin(t=0.1,unit="cm"))    ,axis.title.y= element_text(margin= margin(r=-0,unit="cm"))    ,axis.text.x= element_text(angle=30,hjust=1)  ))to_reciprocal<-function(x) {  ifelse(x>1    , as.character(round(x))    , paste0("1/", round(1/x))  )} ggplot(dat)+  aes(x=n,y=p)+  geom_raster(aes(fill= log(jab_p)),interpolate=TRUE)+# geom_hline(yintercept = p_boundaries, color = "white", alpha = 0.2) +  geom_textcontour(    aes(z=jab_p,label= to_reciprocal(after_stat(level)))    ,color="white"    ,breaks= c(1/30,1/10,1/3,3,10,30,100)  )+  geom_text(    aes(x=n,y=p,label=label,angle=angle)    ,data=evidence_labels    ,color="white"    ,fontface="bold"    ,size=5  )+plot_settings#> Warning: Removed 100 rows containing non-finite values (`stat_textcontour()`).