How do I compute HPDI from anovaBF results?

sherloconan

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

EJ

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.

sherloconan

Thank you, EJ, for an immediate respnse.

Moreover, regarding credible interval, I can sample from anovaBF() object by posterior() 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.

#posterior(): Sample from the posterior distribution of the numerator of a Bayes factor object.
set.seed(123)
chains <- posterior(BF_test, iterations = 10000, progress = FALSE)
summary(chains)

# 1. Empirical mean and standard deviation for each variable,
# plus standard error of the mean:
#   
#                  Mean       SD Naive SE Time-series SE
# mu            12.0140   1.9739 0.019739       0.019739
# Level-Level1  -0.9126   0.2330 0.002330       0.003028
# Level-Level2   0.9126   0.2330 0.002330       0.003028
# Subject-s1    -0.5047   2.0783 0.020783       0.020783
# Subject-s10   -3.4618   2.0686 0.020686       0.020686
# Subject-s2    -4.9396   2.0545 0.020545       0.020545
# Subject-s3     0.4865   2.0692 0.020692       0.020692
# Subject-s4    10.3188   2.0715 0.020715       0.020715
# Subject-s5     4.9143   2.0709 0.020709       0.020709
# Subject-s6     3.9306   2.0735 0.020735       0.020735
# Subject-s7   -10.8390   2.0742 0.020742       0.020742
# Subject-s8     1.4563   2.0734 0.020734       0.020734
# Subject-s9    -1.5008   2.0776 0.020776       0.020776
# sig2           0.9783   0.6581 0.006581       0.018905
# g_Level       12.7965 238.8826 2.388826       2.388826
# g_Subject     53.4141  43.4500 0.434500       0.945808
# 
# 2. Quantiles for each variable:
#   
#                  2.5%      25%      50%     75%    97.5%
# mu             8.0719  10.7764  12.0113 13.2367  16.0214
# Level-Level1  -1.3428  -1.0629  -0.9218 -0.7730  -0.4216
# Level-Level2   0.4216   0.7730   0.9218  1.0629   1.3428
# Subject-s1    -4.6733  -1.8250  -0.5038  0.8030   3.6065
# Subject-s10   -7.6268  -4.7830  -3.4427 -2.1541   0.6773
# Subject-s2    -8.9768  -6.2410  -4.9498 -3.6178  -0.9687
# Subject-s3    -3.7183  -0.8139   0.4960  1.7873   4.6177
# Subject-s4     6.2335   9.0038  10.3119 11.6187  14.4716
# Subject-s5     0.7911   3.6052   4.9123  6.2190   8.9987
# Subject-s6    -0.1594   2.6103   3.9218  5.2237   8.0863
# Subject-s7   -15.0172 -12.1405 -10.8177 -9.5200  -6.7471
# Subject-s8    -2.6950   0.1486   1.4424  2.7655   5.5848
# Subject-s9    -5.7281  -2.7731  -1.4747 -0.2094   2.5858
# sig2           0.3241   0.5555   0.7983  1.1787   2.7453
# g_Level        0.1625   0.7346   1.7170  4.5263  53.7451
# g_Subject      9.1091  25.5156  42.1003 67.6080 167.7259

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".

richarddmorey

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)

sherloconan

Thank you, Dr. Wagenmakers and Dr. Morey, for the responses.