d_waic.R

mod_1 <- brm(value ~ questionnaire * condition, 
             data = idealised_data,
             seed = 42,
             cores = 4)

mod_2 <- brm(value ~ questionnaire + condition, 
             data = idealised_data,
             seed = 42,
             cores = 4)

mod_3 <- brm(value ~ questionnaire, 
             data = idealised_data,
             seed = 42,
             cores = 4)

mod_4 <- brm(value ~ condition, 
             data = idealised_data,
             seed = 42,
             cores = 4)


all_waic <- waic(mod_1, 
                 mod_2, 
                 mod_3, 
                 mod_4)

weight <- function(d_waic, all_d_waic){
  
  num <- exp(-0.5 * d_waic)
  den <- sum(exp(-0.5 * all_d_waic[1]), 
             exp(-0.5 * all_d_waic[2]), 
             exp(-0.5 * all_d_waic[3]), 
             exp(-0.5 * all_d_waic[4]))
  
  return(num / den)
  
}

all_waic <- c(749.7, 755.4, 754.5, 753.6)
d_waic <- all_waic - all_waic[which.min(all_waic)]

# model 1 is superior
weight(d_waic[1], d_waic)
weight(d_waic[2], d_waic)
weight(d_waic[3], d_waic)
weight(d_waic[4], d_waic)

# compare model 1 to model 4 in terms of Kullback–Leibler discrepancy
weight(d_waic[1], d_waic) / weight(d_waic[4], d_waic)

# evidence ratio as the normalized probability that Model 1 is to be 
# preferred over Model 4
weight(d_waic[1], d_waic) / (weight(d_waic[1], d_waic) + 
                               weight(d_waic[4], d_waic))