How do I compute HPDI from anovaBF results?
I am interested in computing highest (posterior) density interval (HPDI / HDI), e.g., in one-way ANOVA within-subjects / repeated-measures design.
I know the scripts to compute a HPDI in Bayesian one-sample t-test, i.e.,
HPDinterval(ttestBF(x = dat_level1, posterior = TRUE, iterations = 100000, progress = FALSE), prob = .95)[1,]
However, HPDinterval() seems not to work on anovaBF() object. The example dataset is attached.
dat_original <- data.frame(
#' Example data (Loftus & Masson, 1994). Long form.
"Subject" = factor(paste("s", rep(1:10,3), sep="")),
"Level" = factor(c(rep("Level1",10),rep("Level2",10),rep("Level3",10))),
"DV" = c(10,6,11,22,16,15,1,12,9,8,
13,8,14,23,18,17,1,15,12,9,
13,8,14,25,20,17,4,17,12,12))
# compare Level1 and Level2 conditions
dat <- subset(dat_original, Level != "Level3")
library(BayesFactor)
set.seed(123)
(BF_test <- anovaBF(DV ~ Level + Subject, data = dat, whichRandom = "Subject",
iterations = 10000, progress = FALSE))
# Bayes factor analysis
# --------------
# [1] Level + Subject : 61.36352 ±1.2%
#
# Against denominator:
# DV ~ Subject
# ---
# Bayes factor type: BFlinearModel, JZS
Comments
I'll attend Richard to your issue (he is busy so may not have time soon -- you could email him and let us know here what he says).
Cheers,
E.J.
Thank you, EJ, for an immediate respnse.
Moreover, regarding credible interval, I can sample from
anovaBF()object byposterior()function. However, using 2.5% and 97.5% quantiles is symmetrical, i.e., equal-tailed interval (ETI). Please, correct me if I misunderstood.By the way, the model is simple as Y_{ij} = \mu_j + b_i + \epsilon_{ij},
where \epsilon_{ij} is i.i.d N(0,\sigma^2) for i = 1,...,N and j=1,...,C; N is the number of subjects, C is the number of levels, Y_{ij} is the response from the i-th subject at the j-th level, \mu_j is the mean of the response at the j-th level, b_i is a mean-zero random effect for the i-th subject.
From the above summary (2. Quantiles for each variable), the 95% credible interval for Level1 and Level2 conditions are respectively [8.0719-1.3428, 16.0214-0.4216] and [8.0719+0.4216, 16.0214+1.3428]. I am not sure whether my interpretation is correct (or perhaps not).
P.S. [Edit-1] grammatical correctness on the original post, "an HPDI".
Looks like you just need to sample from the posterior:
(BF_test <- anovaBF(DV ~ Level + Subject, data = dat, whichRandom = "Subject",
progress = FALSE))
BF_samples = posterior(BF_test, iterations=10000)
HPDinterval(BF_samples)
Thank you, Dr. Wagenmakers and Dr. Morey, for the responses.