Bayes Factor for linear mixed models in R

I am trying to compute the Bayes Factor (BF) for one of the fixed effect with the BayesFactor package in R.

The data has the following structure:

• rating is the dependent variable

• cond is the independent variable with 3 levels ("A", "B", "C")

• C1 is a contrast code derived from cond that opposes "A" (coded -0.50) to "B" and "C" (both coded -0.25)

• C2 is a contrast code derived from cond that opposes "B" (coded -0.50) to "C" (coded +0.5; and "A" is coded 0)

• judge and face are random factors such that face is crossed with judge but nested within cond (and thus also nested within C1 and C2)

DT <- fread("http://matschmitz.github.io/dataLMM.csv")
DT[, judge := factor(judge)]
DT[, face  := factor(face)]

# > DT
#       judge face cond    C1  C2 rating
#    1:    66   13    A -0.50 0.0      1
#    2:    20   13    A -0.50 0.0      4
#    3:    22   13    A -0.50 0.0      7
#    4:    69   13    A -0.50 0.0      1
#    5:     7   13    A -0.50 0.0      3
#   ---                                 
# 4616:    45   62    C  0.25 0.5      2
# 4617:    30   62    C  0.25 0.5      6
# 4618:    18   62    C  0.25 0.5      4
# 4619:    40   62    C  0.25 0.5      3
# 4620:    65   62    C  0.25 0.5      1

Ideally I would like to test the "full" model as in:

library(lmerTest)

lmer(rating ~ C1 + C2 + (1 + C1 + C2|judge) + (1|face), data = DT)

and compute the BF for C1.

I managed to compute the BF for C1 but with random intercepts only:

library(BayesFactor)

BF1 <- lmBF(rating ~ C1 + C2 + judge + face, whichRandom = c("judge", "face"), data = DT)
BF0 <- lmBF(rating ~ C2 + judge + face, whichRandom = c("judge", "face"), data = DT)
BF10 <- BF1 / BF0

# > BF10
# Bayes factor analysis
# --------------
# [1] C1 + C2 + judge + face : 0.4319222 ±15.49%
# 
# Against denominator:
#   rating ~ C2 + judge + face 
# ---
# Bayes factor type: BFlinearModel, JZS

I tried without success this solution to include the random slopes:

BF1 <- lmBF(rating ~ C1 + C2 + judge + face + C1:judge + C2:judge,
            whichRandom = c("judge", "face", "C1:judge", "C2:judge"), data = DT)
# Some NAs were removed from sampling results: 10000 in total.

I would also need to include (if possible) the correlation between the random intercepts and slopes for judge.

Please feel free to use any other package (e.g., rstan, bridgesampling) in your answer.

Some additional questions:

• Do I need to perform any transformation on the BF10, or can I interpret it as it?
• What are the default priors?