Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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,]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.