Include replication Bayes Factor as outcome

[LukasWallrich]

Separating this from the unwieldy #4 which we can hopefully close soon.

Replication Bayes Factors would be a helpful alternative success criterion - introduced here with code available here - but the code is not entirely clear, or neatly packaged, and the copyright is unclear - so this is about more than copy-pasting.

[LukasRoeseler]

This is the R-code Patrick Smela shared:

## Replication Bayes Factor ####

### Prior Setup ####
priorRho <- function(rho, alpha=1) {
  priorDensity <- 2^(1-2*alpha)*(1-rho^2)^(alpha-1)
  logNormalisationConstant <- -lbeta(alpha, alpha)
  result <- exp(logNormalisationConstant)*priorDensity
  return(result)
}

estimationPosteriorU <- function(rho, n, r, alpha=1){
  dataTerm <- (1-rho^2)^((n-1)/2)/((1-rho*r)^((2*n-3)/2))*priorRho(rho, alpha)
  hyperTerm <- Re(hypergeo(1/2, 1/2, (2*n-1)/2, 1/2+1/2*r*rho))
  myResult <- dataTerm*hyperTerm
  return(myResult)
}

estimationPosteriorNormalisationConstant <- function(n, r, alpha=1){
  # The normalisation constant for the replication Bayes factor
  integrand <- function(x){estimationPosteriorU(x, n, r, alpha)}
  myResult <- integrate(integrand, -1, 1)$value
  return(myResult)
}

estimationPosterior <- function(rho, n, r, alpha=1){
  normalisationConstant <- estimationPosteriorNormalisationConstant(n, r, alpha)
  myResult <- 1/normalisationConstant*estimationPosteriorU(rho, n, r, alpha)
  return(myResult)
}

repPrior <- function(rho, nOri, rOri){
  estimationPosterior(rho, n=nOri, r=rOri)
}

### Posterior Setup ####
repPosteriorU <- function(rho, nOri, rOri, nRep, rRep){
  # Unnormalised posterior for the replication Bayes factor
  dataTerm <- (1-rho^2)^((nOri+nRep-2)/2)/((1-rho*rRep)^((2*nRep-3)/2)*(1-rho*rOri)^((2*nOri-3)/2))
  hyperTerm <- Re(hypergeo(1/2, 1/2, (2*nOri-1)/2, 1/2+1/2*rOri*rho))
  myResult <- dataTerm*hyperTerm
  return(myResult)
}

repPosteriorNormalisationConstant <- function(nOri, rOri, nRep, rRep){
  # The normalisation constant for the replication Bayes factor
  integrand <- function(x){repPosteriorU(x, nOri, rOri, nRep, rRep)}
  myResult <- integrate(integrand, -1, 1)$value
  return(myResult)
}

repPosterior <- function(rho, nOri, rOri, nRep, rRep){
  normalisationConstant <- repPosteriorNormalisationConstant(nOri, rOri, nRep, rRep)
  myResult <- 1/normalisationConstant*repPosteriorU(rho, nOri, rOri, nRep, rRep)
  return(myResult)
}

### Bayes Factor Setup ####

repBfR0 <- function(rho=0, nOri, rOri, nRep, rRep){
  myResult <- repPrior(rho, nOri, rOri)/repPosterior(rho, nOri, rOri, nRep, rRep)
  return(myResult)
}

### Bayes Factor Calculation ####

BFR0 = NULL # Create Dummy to store output

# Unfortunately the Bayes factor calculation contains a bug that occurs when certain parameters become too large. Here is a workaround
N_workaround <- df$N_Repl                # Copy replication sample sizes into a dummy
N_workaround[N_workaround > 1200] = 1200 # Set all replication sample sizes larger 1200 to 1200
r_workaround <- df$r_Repl                # Copy replication correlation coefficients into a dummy
r_workaround[r_workaround > .7] = .7     # Set all replication correlation coefficients larger .7 to .7

# Run Loop Analysis over each study
for (i in 1:length(df$N_Original)) {
  
  BFR0[i] <- repBfR0(nOri = df$N_Original[i], # Sample Size of the original study
                     rOri = df$r_Original[i], # Effect Size of the original study as correlation
                     nRep = N_workaround[i],  # Sample Size of the replication study
                     rRep = r_workaround[i])  # Effect Size of the replication study as correlation
  
}

[LukasRoeseler]

Dennis Kondzic talked about measuring replication success at DGPs 2024 using difference and equivalence tests that we may want to incorporate into a more exhaustive discussion of success criteria:

[LukasWallrich]

I wonder how many of these are the same thing with different directions (i.e. difference vs equivalence) - but I like the idea of clearly outlining what questions one might want to answer, and what test is appropriate.

…

On Wed, 18 Sept 2024, 09:15 Lukas Röseler, ***@***.***> wrote: Dennis Kondzic talked about measuring replication success at DGPs 2024 using difference and equivalence tests that we may want to incorporate into a more exhaustive discussion of success criteria: Screenshot.2024-09-18.100634.png (view on web) <https://github.com/user-attachments/assets/1de39fca-1dfe-4e1b-9521-a2fbdacdc15a> Screenshot.2024-09-18.100726.png (view on web) <https://github.com/user-attachments/assets/eb3d1393-8870-47fc-be0d-bfc28e44be9a> image.png (view on web) <https://github.com/user-attachments/assets/3767d155-1535-4a33-9a8b-fe4bf5ad7b99> — Reply to this email directly, view it on GitHub <#50 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AOK6NGN4VMOFJYDL4ELWWDTZXEZBLAVCNFSM6AAAAABMORN3ISVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNJXHAYDGMJQHA> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[LukasRoeseler]

This may help us including bayesian factors: https://cran.r-project.org/web/packages/ReplicationSuccess/index.html

Edit: Another bayesian criterion: Verhagen, J., & Wagenmakers, E.-J. (2014). Bayesian tests to quantify the result of a replication attempt. Journal of Experimental Psychology: General, 143(4), 1457–1475. https://doi.org/10.1037/a0036731

[LukasRoeseler]

This may help us aggregating bayes factors: https://doi.org/10.1002/jrsm.1765