Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Supported by

Model averaged Q-Q plots in R

Hi all,

I would like to reproduce the model averaged Q-Q plots that can be plotted by selecting the "Q-Q plots of residuals" option under the Bayesian (repeated measures) ANOVA in JASP. I'm using the BayesFactor package of Richard Morey in R.

Is there a way to obtain the code used to generate these plots, or alternatively an explanation on how to obtain the residuals?

Many thanks,

Myrthel

Comments

  • I'll ask Don van den Bergh, who programmed this

  • edited February 19

    Unfortunately this is not trivial to implement.

    The relevant functions are .BANOVAsamplePosteriors (https://github.com/jasp-stats/jasp-desktop/blob/development/JASP-Engine/JASP/R/commonAnovaBayesian.R#L1506) and .BANOVAgetSMIResidRsq (https://github.com/jasp-stats/jasp-desktop/blob/development/JASP-Engine/JASP/R/commonAnovaBayesian.R#L2003). However, to use those you need the complete JASP environment. In addition, we do some tricks to avoid storing all posterior samples for all models in memory at once, which not be a problem in your situation.


    A small example in R with two models:


    library(BayesFactor)

    data("puzzles")


    a1 <- anovaBF(RT ~ shape + color, data = puzzles, progress=FALSE)


    samples_1 <- posterior(a1, 1, iterations = 1e3)

    samples_2 <- posterior(a1, 2, iterations = 1e3)


    samples_1 <- as.matrix(samples_1)[, 1:3]

    samples_2 <- as.matrix(samples_2)[, 1:3]


    datOneHot_1 <- model.matrix(RT ~ shape, data = puzzles[,c("RT", "shape")], contrasts.arg = lapply(puzzles[,"shape", drop = FALSE], contrasts, contrasts = FALSE))

    datOneHot_2 <- model.matrix(RT ~ color, data = puzzles[,c("RT", "color")], contrasts.arg = lapply(puzzles[,"color", drop = FALSE], contrasts, contrasts = FALSE))


    preds_1 <- tcrossprod(datOneHot_1, samples_1)

    preds_2 <- tcrossprod(datOneHot_2, samples_2)


    resid_1 <- puzzles[["RT"]] - preds_1

    resid_2 <- puzzles[["RT"]] - preds_2


    gives you the residuals for two models. Next you need to weigh these according to the posterior model probabilities, after which you get the model averaged distribution over residuals. Using that distributions you can then compute things like credible intervals or make plots.


    If you have any more questions, please let me know!

Sign In or Register to comment.