Probability-of-Bivariate-Superiority (PBS) Estimates (R Code)

Li, J. C.-H., & Waisman, R. (2019). Probability of bivariate superiority: A non-parametric common-language statistic for detecting bivariate relationships. Behavior Research Methods, 51, 258-279. https://doi.org/10.3758/s13428-018-1089-5

Description: Obtaining Estimates based on Probability-of-Bivariate-Superiority

Copy and paste the R code below, or get access to the R file via: https://www.dropbox.com/scl/fi/qfefr09oxyvqu76v1kigo/pbs_function.R?rlkey=j2mo4kgjmkwxaz3vru61y5e6x&st=45hwl7rz&dl=0

• #’ A pbs Function
• #’
• #’ This function allows you to obtain estimates based on Probability-of-Bivariate-Superiority
• #’ @keywords pbs
• #’ @export
• #’ @examples
• #’ The following example shows the details of how to obtain the estimates based on probability-of-bivariate-superioirty.
• #’
• #’ This example is based on:
• #’ Li, J. C.-H., & Waisman, R. (2019). Probability of bivariate superiority: A non-parametric 
• #’ common-language statistic for detecting bivariate relationships. Behavior Research 
• #’ Methods, 51, 258-279. https://doi.org/10.3758/s13428-018-1089-5
• #’
• #’ First, you can save the following scores as a 100 by 2 data matrix in R or RStudio
• #’X <-  c(1.560, .032, 1.648, -1.040, -.067, -1.231, 1.669, .803, -.005, -1.476, -.383,
• #’-.473, 1.620, -1.082, -1.088, -1.515, -1.717, -1.642, .651, -.587, -.608, -.238, .385,
• #’.626, 1.625, .585, .807, -.811, 1.653, 1.106, -.629, -1.651, 1.104, 1.074, -.961, .805,
• #’-1.086, 1.124, 1.705, .810, 1.338, .283, -.480, .527, .930, -.158, -.817, 1.645, 1.612,
• #’-1.603, .853, -.065, .605, .742, 1.102, 1.631, -1.631, -.678, .740, .602, .149, -1.097,
• #’-.907, -1.536, -.263, -1.632, -.844, 1.330,        -.493, -.983, -.068, 1.130, -.589,
• #’1.137, -1.267, -.842, -.069, 1.568, 1.691, .645, -.378, .526, -.740, -1.204, -.777,
• #’-1.572, -1.392, -1.707, -1.228, -.720, -1.067, -.153, .260, -1.173, -.636, -.177, .457, .222, -.574, -.966)
• #’
• #’Y  <- c(1.086, -1.685, -.012, 1.624, .162, -1.281, -1.203, -1.064, 1.522, -.596, .568,
• #’-.256, 1.239, -.301, 1.040, -1.387, .828, -1.237, .605, 1.449, .949, 1.268, 1.500, 1.591,
• #’-1.579, 1.065, -1.458, .702, 1.380, -.387, -.060, 1.637, .642, 1.156, 1.372, -.579, -1.571,
• #’-.343, .520, -.114, .192, -1.352, .745, -1.338, .672, -1.632, -.866, .674, -1.548, 1.182, .401,
• #’-1.109, .084, .400, .660, 1.709, -1.605, -.291, .423, 1.208, .132, -.049, -.921, -1.126, -.382,
• #’.366, 1.024, -1.428, -.114, -1.672, -.177, 1.103, -.063, .610, -.886, 1.138, -1.325, .530, -.915,
• #’-1.455, .424, .688, .889, -.937, -1.514,         -.343, -.612, -1.062, -.137, -1.244, .542, -.191,
• #’-1.191, 1.038, -.973, 1.140, 1.493, .751, 1.402, -.074)
• #’data <- matrix(c(X,Y),100,2)
• #’
• #’Second, you can execute the following function
• #’pbs(data,1000,.95, 1234, 4)
• #’
• #’Third, you can obtain the results below
• #’                         Results
• #’Bp                        0.6100
• #’BSI lower limit           0.5031
• #’BSI upper limit           0.7169
• #’BPI lower limit           0.4800
• #’BPI upper limit           0.7000
• #’BCaI lower limit          0.5100
• #’BCaI upper limit          0.7200
• #’Bootstrap Percentage %        95
• #’No. of bootstrap samples    1000
• #’Seed no.                    1234
• pbs <- function(data, bootno = 1000, percent = .95, seed.no = 1234, dec = 4){
•   require(MASS)
•   require(boot)
•   set.seed(seed.no)
•   pb<-function(data){
•     n=length(data[,1])
•     c=n
•     tcount2=0
•     meanx=mean(data[,1])
•     meany=mean(data[,2])
•     for(i in 1:n){
•       diff2=(data[i,1]-meanx)*(data[i,2]-meany)
•       if(diff2>0){tcount2=tcount2+1}
•       if(diff2==0){tcount2=tcount2+0.5}}
•     return(tcount2/(c))}
•   spec <- function(x, dec) format(round(x, dec), nsmall=dec)
•   funbp <- function(d, i){
•     d2 <- d[i,]
•     return(pb(d2))}
•   bootr <- boot(data, funbp, R = bootno)
•   bci <- boot.ci(bootr, conf = percent, type = c(“perc”,”bca”))
•   sebsi <- sd(bootr[[2]])
•   bsill <- qnorm((1-percent)/2)*sebsi+bci[[2]]
•   bsiul <- qnorm((1-percent)/2+percent)*sebsi+bci[[2]]
•   re <- array(0, dim = c(10,1))
•   rownames(re) <- c(“Bp”,”BSI lower limit”,”BSI upper limit”,”BPI lower limit”,
•                     “BPI upper limit”, “BCaI lower limit”, “BCaI upper limit”, “Bootstrap Percentage %”,
•                     “No. of bootstrap samples”, “Seed no.”)
•   colnames(re) <- “Results”
•   re[1,1] <- spec(bci[[2]],dec)
•   re[2,1] <- spec(bsill,dec)
•   re[3,1] <- spec(bsiul, dec)
•   re[4,1] <- spec(bci[[4]][4], dec)
•   re[5,1] <- spec(bci[[4]][5], dec)
•   re[6,1] <- spec(bci[[5]][4], dec)
•   re[7,1] <- spec(bci[[5]][5], dec)
•   re[8,1] <-  percent*100
•   re[9,1] <- bootno
•   re[10,1] <- seed.no
•   re <- as.data.frame(re)
•   return(re)}