Random slopes in lmBF

[jen265]

Hello, I am trying to use the lmBF function with random slopes, but it always results in the following error:

In doNwaySampling(method, y, X, rscale, iterations, gMap, incCont, : Some NAs were removed from sampling results: 10000 in total.

The dataset does not have missing values. The model includes Subject as the grouping variable, IV1 as the level-1 predictor with a random slope, and IV2 to IV6 as level-2 predictors, and is specified as follows:

full_BF = lmBF(DV ~ IV1 IV2 + IV1 IV3 + IV1 IV4 + IV1 IV5 + IV1 * IV6 + Subject + IV1:Subject, data = data, whichRandom = c('Subject','IV1:Subject'))

Any help to solve the error is very much appreciated.

[richarddmorey]

Can you provide a reproducible example?

[jen265]

Thank you for your swift reply. Here is an example:

install.packages(c("data.table","BayesFactor"))

library(data.table) library(BayesFactor) data <- fread("http://eng1.psy.de/data.csv") data$Subject2 <-as.factor(data$Subject) data$Subject <- NULL names(data)[names(data) == "Subject2"] <- "Subject" full_BF = lmBF(DV ~ IV1 IV2 + IV1 IV3 + IV1 IV4 + IV1 IV5 + IV1 * IV6 + Subject + IV1:Subject, data = data, whichRandom = c('Subject','IV1:Subject')) full_BF

[richarddmorey]

I've confirmed the issue but not sure what's causing it. I'll get back to you when I've explored it.

[richarddmorey]

How were these data constructed? I notice that the mean IV1 for each subject is (computationally) 0, so I suspect this issue is arising because there are too many parameters in the model. This could be caused by implicitly eliminating them through centering, etc.

[jen265]

The mean IV1 for each subject is 0 because IV1 is group-mean centered. Moreover, all variables are z-standardized. As we are interested in the cross-level interactions, group-mean centering of the level-1 predictor seems necessary. When we analyze the data with the lmer function, the model is computed without errors.

EDIT: The error also occurs when I use the raw data (no variables centered or z-standardized) or the following simplified model: full_BF = lmBF(DV ~ IV1 * IV2 + Subject + IV1:Subject, data = data, whichRandom = c('Subject','IV1:Subject'))

[jen265]

Has someone already succeeded in using the lmBF funtion with random slopes? I tried it with different datasets and it always resulted in the known error.

[jen265]

A manuscript on the BayesRS package (https://psyarxiv.com/4xqvr/) says: „The BayesFactor package (Morey & Rouder, 2014), which computes BFs for mixed models, allows users to specify random intercepts and random slopes for categorical predictors (i.e., variables on a nominal scale), as are often used in experimental research. However, at present, there is no statistical software that allows users to readily compute BFs for models that include random slopes on continuous predictors“ (p. 7).

Can the error described in my first post be attributed to the fact that my level-1 predictor is continuous? When I recode my continuous level-1 predictor into a categorial predictor (only 0 and 1 as values), the error does not occur. So would it be correct to write in a manuscript that the BayesFactor packgage has not allowed users to specify random slopes for continuous predictors yet?

[m-Py]

Sorry for highjacking this thread. After reading a new preprint ("Bayes Factors for Mixed Models": https://psyarxiv.com/y65h8) that uses BayesFactor to compute Bayes factors for mixed models (including the random slope of a condition over participants), I also tried out BayesFactor for this purpose. I ran into the same problem as @jen265 with a data set of my own. It is also possible reproduce the behaviour using a sample data set from the package lme4, which may help tackle the problem, if you find any time to check it out:

BayesFactor::lmBF( 
  Reaction ~ Days * Subject,  
  data = lme4::sleepstudy,  
  whichRandom = "Subject"
)

Anyway, thanks for developing BayesFactor, it is greatly appreciated!

EDIT: I can also confirm that the issue does not occur when the predictor is a factor.

[richarddmorey]

Thanks for sharing the example with the simpler data set. My guess is that this is a model constraint issue. If I do:

# Use development version to fix weird character issue
devtools::install_github("richarddmorey/BayesFactor", subdir = "pkg/BayesFactor")

library(BayesFactor)

ss = lme4::sleepstudy
ss$Days = ss$Days - 5 # help estimation of intercept; reduces sensitivity on prior so we can identify issue

bf = BayesFactor::lmBF( 
  Reaction ~ Subject + Days:Subject,  
  data = ss,
  whichRandom = "Subject"
) 

This model gives each Days:Subject interaction its own slope. with a random intercept. It only fails when I add the Days "average" slope, which gives me too many parameters. This has to do with the way these models (and the interactions) are defined (see Rouder et al, 2012)

I can estimate the equivalent model in lme4 (no correlation between slope and intercept, random slope and intercept) and the estimates are essentially identical.

ch = posterior(bf, iterations = 10000)
bf_intercept = colMeans(ch[,1+1:18] + ch[,1])
bf_slopes = colMeans(ch[, 1 + 18 + 1:18])

lme = lme4::lmer(Reaction ~ (-1+Days|Subject) + (1|Subject), data = ss)

plot(bf_intercept, coef(lme)$Subject$`(Intercept)`)
abline(0,1)
plot(bf_slopes, coef(lme)$Subject$Days)
abline(0,1)

If you add in the fixed effects of Days, this will change where the model does pooling. To replicate that, we'd need to add constraint on the random slopes to "pop" the fixed effect out. We do this by treating subjects as "fixed" so that we get the right number of slope parameters.

(here I've added in a separate "Subject0" variable to retain the random intercept parametrization)

ss$Subject0 = ss$Subject

bf2 = BayesFactor::lmBF( 
  Reaction ~ Days + Subject0 + Days:Subject,  
  data = ss,
  whichRandom = "Subject0"
) 

ch2 = posterior(bf2, iterations = 10000)
bf_intercept2 = colMeans(ch2[,2+1:18] + ch2[,1])
bf_slopes2 = colMeans(ch2[, 2 + 18 + 1:18] + ch2[,2])

lme2 = lme4::lmer(Reaction ~ Days + (-1+Days|Subject) + (1|Subject), data = ss)

plot(bf_intercept2, coef(lme2)$Subject$`(Intercept)`)
abline(0,1)
plot(bf_slopes2, coef(lme2)$Subject$Days)
abline(0,1)
abline(h = lme4::fixef(lme2)[2])
abline(v = mean(ch2[,"Days-Days"]))

The estimates are very close, with different amounts of pooling due to different priors.

[m-Py]

Thanks for your elaborate response, very educational. So let's assume I would like to compare a model that does not assume a random slope with a model that does, the following might work?

library(BayesFactor)
library(lme4)
ss = sleepstudy
ss$Subject0 = ss$Subject

bf1 = BayesFactor::lmBF( 
  Reaction ~ Days + Subject,
  data = ss,
  whichRandom = "Subject"
) 

bf2 = BayesFactor::lmBF( 
  Reaction ~ Days + Subject + Days:Subject0,
  data = ss,
  whichRandom = "Subject"
) 

bf2 / bf1
#> Bayes factor analysis
#> --------------
#> [1] Days + Subject + Days:Subject0 : 561.3099 ±1.58%
#> 
#> Against denominator:
#>   Reaction ~ Days + Subject 
#> ---
#> Bayes factor type: BFlinearModel, JZS

Here, the model assuming the random slope is favored by a factor of about 560. This is somewhat consistent with a different kind of Bayes factor that Paul Bürkner computed using brms for the same model comparison (see https://osf.io/m493h/), even though that Bayes factor was larger by an order of magnitude, which I assume may be due to different priors.

[richarddmorey]

That looks right. I'd be careful interpreting a high BF from brms, because I suspect that the bridge sampler (which brms uses) will not be as accurate as the custom-tailored importance sampler that BayesFactor uses. The bridge sampler is more flexible but will have accuracy difficulty in the tails (though, to be fair, that's where it typically doesn't matter). I imagine that the differences between brms and BayesFactor are a combination of sampler and prior differences.

[m-Py]

Thanks a lot for your input, Richard!

[amh571]

Hi! Sorry also for hi-jacking this thread and for (potentially) a basic question. I was wondering what would be the difference between including random slopes with the notation:

e.g. RT ~ condition + subject + 1:subject, whichRandom = "subject"

and

e.g. RT ~ condition + subject + condition:subject, whichRandom = "subject" ?

Would this have to do with the assumed hierarchical structure of the random effects (crossed or nested)? So far all the examples I have seen use the latter. Thank you!

[richarddmorey]

The notation

RT ~ condition + subject + 1:subject, whichRandom = "subject"

isn't supported - by including + subject in the formula you've introduced the random intercent of subject (assuming whichRandom includes subject as well). Your second formula is correct, if you want to include the random intercept and slope.

[amh571]

Thank you for the quick reply and for the explanation! Could you please detail more on what do you mean by 'not supported'? Would it mean that the first formula essentially includes only intercepts, and therefore there would be no difference between:

RT ~ condition + subject + 1:subject, whichRandom = "subject"

and

RT ~ condition + subject, whichRandom = "subject"?

NB: just to clarify, I want to include both random intercepts and slopes, and the data is crossed rather than nested.

[richarddmorey]

Interactions are built by multiplying terms together, so 1:subject will create terms that are all 1*subject, which is the same as subject. By including both subject and 1:subject you've included two identical terms.

By including subject in the formula you've told BayesFactor it is an additive term, then by including subject in whichRandom you've indicated it is random. A random additive term is a random intercept.

It becomes a random slope when it interacts with other terms.

[Catsine]

Hello,

apologies for hacking this thread as well. I am trying to create a model in which I have a random slope for two main effects and their interaction, while the random factor is the participants in the experiment. My model looks like this: BF_full = lmBF(RTs~condition*task+X1+X2+X3+ID+condition*task:ID, whichRandom=c('ID'),data=data) which I compare to: BF_noint = lmBF(RTs~condition+task+X1+X2+X3+ID+condition*task:ID, whichRandom=c('ID'),data=data)

I thought that this should be the right way to write the model, but when I do compare it to the results provided by a lme model comparison (Log-lik ratio test), the two give me point me towards opposite directions. While for the lme model the difference is significant and I should include the interaction term, from the bayesian analysis, I get a bayes factor of 0.20, which slight evidence for the null hypothesis.

Is the way I formulated the model with the lmBF function correct or is there a mistake? I tried to look for examples online dealing with an interaction term in the random slope, but I could not find any.

Many thanks!

[richarddmorey]

Hi @Catsine, there's nothing wrong with your notation per se but the main question is whether that's what you want to test (and what about the data might be driving the difference, if it is a difference). Plots, etc would help but aren't really appropriate here. You might try the forum here: https://forum.cogsci.nl/categories/jasp-bayesfactor

[Catsine]

Many thanks @richarddmorey, I think in the end that the results make sense now. I was running few iterations and the variance of the output was very high, but now that I increased them I get results that are in line with the frequentist analysis. Thanks again!