This is my attempt at an `emmeans` integration for `BayesFactor` R package

BFBayesFactor_2_emmeans.R

#' @export
recover_data.BFBayesFactor <- function (object,index = 1, ...) {
  object <- object[index]
  
  data <- object@data
  formula <- as.formula(object@numerator[[1]]@longName)
  trms <- delete.response(terms(eval(formula, parent.frame())))
  
  # Remove random effects
  random_ <- object@numerator[[1]]@dataTypes=='random'
  if (any(random_)) {
    data <- data[,colnames(data)!=names(random_)[random_]]
    formula <- as.formula(paste0('~. -',names(object@numerator[[1]]@dataTypes)[random_],':.'))
    trms <- update(trms,formula)
  }
  
  cl <- call("mcmc.proxy", formula = formula, data = quote(data))
  emmeans::recover_data(cl, trms, NULL, data, ...)
}

#' @export
#' @import BayesFactor
emm_basis.BFBayesFactor <- function (object, trms, xlev, grid,
                                     iterations = 10000,seed = NULL,index = 1,
                                     est_intercept,...){
  set.seed(seed)
  object <- object[index]
  
  # model.matrix
  m <- model.frame(trms, grid, na.action = na.pass, xlev = xlev)
  factors_ <- sapply(m, is.factor)
  if (any(factors_)){
    m[factors_] <- lapply(m[factors_], make_full_dummy)
  }
  X <- model.matrix(trms, m)
  
  # Re-compute intercept from mu
  if (missing(est_intercept)) {
    post_samples <- as.data.frame(posterior(object,iterations = iterations,progress = FALSE))
    
    post_samples <- trim_rename_sample_frame(post_samples,object, trms, xlev, X)
    
    A <- suppressMessages(emmeans::emmeans(object,~1,est_intercept = post_samples))
    A <- A@linfct %*% t(A@post.beta)
    post_samples[,1] <- post_samples[,1] - A
  } else {
    post_samples = est_intercept;
    post_samples[,1] <- 0
  }
  
  # Bhat and V
  bhat = apply(post_samples, 2, median)
  V = cov(post_samples)
  
  misc = list(model = object)
  
  list(X = X,
       bhat = bhat,nbasis = matrix(NA),V = V,
       dffun = function(k,dfargs) Inf,dfargs = list(),
       misc = misc,post.beta = post_samples)
}

make_full_dummy <- function(fct, is.random = FALSE){
  CC <- stats:::.Diag(levels(fct),FALSE)
  if (is.random){
    CC[] <- 0
  }
  contrasts(fct,how.many = length(unique(fct))) <- CC
  fct
}

trim_rename_sample_frame <- function(sample_frame,object, model_trms, xlev, model_matrix){
  data <- object@data
  formula <- as.formula(object@numerator[[1]]@longName)
  # formula <- flip_interactions(formula)
  
  trms <- delete.response(terms(eval(formula, parent.frame())))
  m <- model.frame(trms,data = data, na.action = na.pass, xlev = xlev)
  factors_ <- sapply(m, is.factor)
  
  if (any(factors_)){
    m[factors_] <- lapply(m[factors_], make_full_dummy)
  }
  X <- model.matrix(trms, m)
  
  sample_frame <- sample_frame[, seq_len(ncol(X)), drop = FALSE]
  
  factor_grid <- attr(trms,"factors")
  are_ints <- factor_grid %>%
    as.data.frame() %>%
    map_lgl(~(sum(.x)>1)) %>%
    which()
  
  sample_frame_old <- sample_frame
  for (i in seq_along(are_ints)) {
    sub_fcts <- factor_grid[,are_ints[i]]
    sub_fcts <- sub_fcts[sub_fcts==1]
    
    sub_fcts_code <- map_dbl(attr(X,"contrasts"),ncol)
    sub_fcts_code <- ifelse(names(sub_fcts) %in% names(sub_fcts_code),
                            sub_fcts_code[names(sub_fcts)],
                            1)
    names(sub_fcts_code) <- names(sub_fcts)
    
    BF_order <- map(sub_fcts_code,seq_len) %>%
      expand.grid() %>%
      as.data.frame() %>%
      rev() %>%
      reduce(interaction) %>%
      as.numeric()
    
    sample_frame[,attr(X,"assign")==are_ints[i]] <-
      sample_frame[,attr(X,"assign")==are_ints[i],drop = FALSE][,BF_order,drop = FALSE]
  }
  
  
  colnames(sample_frame) <- colnames(X)
  sample_frame <- sample_frame[, colnames(sample_frame) %in% colnames(model_matrix), drop = FALSE]
  
  return(as.matrix(sample_frame))
}