gistfile1.R

library(magrittr)
library(dplyr)
library(BayesFactor)

theseData <- data.frame(
  estimate = c(104.452368548111,117.413178585352,78.7223741725136,76.4106242567969,82.7631936151179,56.9203980838411,74.4102158311052,87.3611057090617,60.9334108059899,119.71391794231,106.264246216,108.228558111754,96.7640506225208,60.7934578606003,42.4173983045924,129.926303462432,13.5213061083478,73.0941378407556,40.4408225175657,86.1817712338958,70.581097150017,48.8449615326597,90.6359161067542,42.2265556453093,98.7564777546867,96.176302161205,62.3612161520409,113.76179390849,46.7560720387586,52.2189402515585,73.6369606772811,106.759109436828,51.0453162114017,94.9706496538741,102.540298026013,60.379138518731,89.6665186074092,93.1648279908775,131.023471716353,132.861306717587,147.252946161977,131.367061791528,146.642024628095,70.8522322296042,92.8410392664767,81.6321679489163,115.498166140208,103.024491359749,104.771202195881,100.629398797581,70.7405794838797,92.2778430688292,76.7296920976546,72.1162388510765,72.4319815851055,146.751555086237,78.0052719819749,216.243937309623,40.8595320101227,80.45287267714,53.2766678111939,113.757326522442,117.107022430137,157.42584783562,131.166605792994,92.5577507743433,130.276701828596,150.031507999395,88.3446557287025,86.1433246394917,107.548712533477,112.417949522607,71.8418647319083,121.176336589237,82.5376641187226,181.359415772665,69.0825006890801,99.1453523196733,97.6061856581658,76.3026377563864,162.613906450788,101.52986940189,214.684409090797,65.2239980303493),
  eccentricity = c(3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5,11.5),
  crowded = c("Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Bouma", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard", "Standard"), 
  ID = c("A" ,"B" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"J" ,"K" ,"L" ,"M" ,"N" ,"O" ,"A" ,"B" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"J" ,"K" ,"L" ,"M" ,"N" ,"O" ,"A" ,"B" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"J" ,"K" ,"L" ,"M" ,"N" ,"O" ,"A" ,"B" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"J" ,"K" ,"L" ,"M" ,"N" ,"O" ,"A" ,"B" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"L" ,"N" ,"O" ,"A" ,"C" ,"D" ,"E" ,"F" ,"G" ,"H" ,"I" ,"K" ,"L" ,"M" ,"N")
)


inclusionBF <- function(priorProbs, variable){
  
  ###https://www.cogsci.nl/blog/interpreting-bayesian-repeated-measures-in-jasp###
  
  
  if(typeof(priorProbs) == 'S4') priorProbs <- as.vector(priorProbs)
  
  
  theseNames <- names(priorProbs)
  nProbs <- 1:length(priorProbs)
  variableMatches <- grep(variable, theseNames)
  
  
  if(grepl(':', variable)){
    subordinateVariables <- variable %>% strsplit(':') %>% unlist()
    
    thisRegex <- paste0(subordinateVariables,collapse = '.*\\+.*')
    
    subordinateEffects <- grep(thisRegex, theseNames, perl = T)
    
    subordinateEffects <- subordinateEffects[!subordinateEffects %in% variableMatches]
    
    
    sum(priorProbs[variableMatches])/sum(priorProbs[subordinateEffects])
  } else {
    interactionMatches <- grep(paste0(variable,'(?=:)|(?<=:)',variable), theseNames, perl = T)
    
    variableMainEffects <- variableMatches[!variableMatches %in% interactionMatches]
    
    
    otherMainEffects <- nProbs[!nProbs %in% c(variableMainEffects,interactionMatches)]
    
    
    sum(priorProbs[variableMainEffects])/sum(priorProbs[otherMainEffects])
  }
}




thisTest <-  generalTestBF(estimate ~ eccentricity * crowded * ID, 
                           data=theseData,
                           whichRandom = c('ID', 'eccentricity:ID', 'crowded:ID', 'eccentricity:crowded:ID'), 
                           progress = TRUE
                        ) 


thisTest

thesePriorProbs <- thisTest %>% newPriorOdds() %>% `*`(thisTest) %>% as.BFprobability() %>% as.vector() #extract P(M|Data) from BF object

theseInclusionBFs <- data.frame(
  variable = c('crowded','eccentricity','eccentricity:crowded','ID', 'eccentricity:ID', 'crowded:ID', 'eccentricity:crowded:ID'),
  BF = -999,
  inclusionEvidence = character(1),
  stringsAsFactors = F
)

for(thisVariable in theseInclusionBFs$variable){ #calculate inclusion BFs for main effects and interactions
  thisInclusionBF <- inclusionBF(thesePriorProbs, thisVariable)
  theseInclusionBFs %<>% mutate(BF = replace(BF, variable == thisVariable, thisInclusionBF))
}

theseInclusionBFs %<>% 
  mutate(inclusionEvidence = replace(inclusionEvidence, BF>3, 'Include')) %>%
  mutate(inclusionEvidence = replace(inclusionEvidence, BF<1/3, 'Exclude')) %>%
  mutate(inclusionEvidence = replace(inclusionEvidence, BF > 1/3 & BF < 3, 'Equivocal'))

thisPosterior <- posterior(thisTest[6], iterations = 10000, progress = FALSE)
thisPosteriorSummary <- summary(thisPosterior)


thisPosteriorSummary$statistics[1:2,] 

thisPosteriorSummary$quantiles[1:2,]