Sampling from Priors

[mattansb]

Hi Richard,

Is it currently possible to get samples for the prior distributions?
(What I would really like to do it compute a BF for a specific "contrast" using the Savage-Dickey density ratio method. I've read you post of equality constraints, but found the application very difficult when trying to break down an interaction by comparing two simple effects).

Thanks!

[richarddmorey]

Not currently, but given that they are just multivariate-Cauchy-distributed (as detailed in Rouder et al 2012) you could just sample directly and then multiply by the matrix of contrast codes.

[mattansb]

Using the g parameter as the scaling of the Cauchy distribution?

For example:

library(BayesFactor)

data(puzzles)
result <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                 progress=FALSE)

chains <- as.data.frame(posterior(result,index = 4, iterations = 10000))

prior_samples_shape <- rcauchy(10000,scale = mean(chains$g_shape))
post_samples_shape_round <- chains$`shape-round`

# test if round parameter is different than 0?
SavageDickeyBF <- approxfun(density(post_samples_shape_round))(0) /
  approxfun(density(prior_samples_shape))(0)

Thanks!

[mattansb]

Sorry, should the scale of the cauchy be sig2? Or both sig2 and g_*?

[richarddmorey]

r is the scale of the Cauchy distribution. The resulting samples will be samples from the standardised contrast effect size (that is, the effect divided by the error standard deviation). But if you're using Savage-Dickey, you don't need the samples themselves; you only need the density at 0, which you can calculate analytically.

[mattansb]

Perhaps I'm missing something - the posterior samples don't seem to be standerdized - that is, they are affected by the scale of the DV:

library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************

data(puzzles)

result <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                  progress=FALSE)
chains <- as.data.frame(posterior(result,index = 4, iterations = 10000))

puzzles$RT <- puzzles$RT*100
result2 <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                   progress=FALSE)

chains2 <- as.data.frame(posterior(result2,index = 4, iterations = 10000))

plot(density(chains$`shape-round` - chains$`shape-square`))

plot(density(chains2$`shape-round` - chains2$`shape-square`))

c(desity_at_0_orig  = approxfun(density(chains$`shape-round` - chains$`shape-square`))(0),
  desity_at_0_trans = approxfun(density(chains2$`shape-round` - chains2$`shape-square`))(0))
#>  desity_at_0_orig desity_at_0_trans 
#>      0.0817700892      0.0007774765

^{Created on 2019-02-25 by the reprex package (v0.2.1)}

[mattansb]

But if you're using Savage-Dickey, you don't need the samples themselves; you only need the density at 0, which you can calculate analytically.

Ah, yes - good point !

[mattansb]

Sorry - I completely misunderstood which part of my message you were replying to!

If I understand correctly:

• The density at 0 of the prior Cauchy distribution should be with dcauchy(0,scale = rscaleFixed)
• And the density at 0 of the posterior should be estimated by standardizing the contrast by dividing by sqrt(sig2).

Like so:

library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************

data(puzzles)
result <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                  progress=FALSE)

chains <- as.data.frame(posterior(result,index = 4, iterations = 10000))
post_samples_shape_round <- with(chains,(`shape-square` - `shape-round`)/sqrt(sig2))

# test if contrast is different than 0? (multiply scale by 2 for the contrast)
SavageDickeyBF <- dcauchy(0,scale = .5*2) / approxfun(density(post_samples_shape_round))(0)
SavageDickeyBF
#> [1] 2.908753

^{Created on 2019-02-25 by the reprex package (v0.2.1)}

Thanks!

[richarddmorey]

Yes, that's one way to do it. Here's my example:

set.seed(1560)

N = 100
k = 2
dat = data.frame(grp = factor(rep(1:k, N)), y = rnorm(N*k))

# for the example
rscale = 1

bf = BayesFactor::lmBF(y ~ grp, data = dat, rscaleFixed = rscale)
chains = BayesFactor::posterior(bf, unreduce = FALSE, iterations = 10000)

# The Rouder et al (2012) priors are based on constrains formed by 
# projecting k - 1 means into k dimensions. In reverse, these are essentially 
# contrasts. "unreduce" above gives us the unprojected posterior, which we need.
# with k=2, this will be a scaled difference between the two conditions
eff_size = chains[,"grp_redu_1"]
# See Rouder et al for more on how this is done.

# logspline estimate
ls_obj = logspline::logspline(eff_size)

curve(logspline::dlogspline(x, ls_obj), min(eff_size), max(eff_size), col = "blue")  
curve(dcauchy(x, scale = rscale), col = "red", add = TRUE)
abline(v = 0, lty = 2)

# Add prior
points(0, dcauchy(0, scale = rscale), pch = 19, col = "red")

# estimate of density at 0, for Savage-Dickey
est0 = logspline::dlogspline(0, ls_obj)
# add point
points(0, est0, col = "blue", pch = 19)

# Savage-Dickey estimate of BF
dcauchy(0, scale = rscale) / est0
# Gaussian quadrature estimate of BF
as.vector(bf)[1]

[richarddmorey]

(Note that in my example I start with the standardised posterior samples, and your example starts with the unstandardised ones and "reduces"/standardises them )

[mattansb]

Thanks Richard!

[mattansb]

If anyone ever ends up here - I've implemented this method in BFEffects::inferBF.

Thanks again!

[mattansb]

One more question regarding priors:
When modeling an interaction between a factor and a covariate (using generalTestBF), is the scale of the prior of the interaction term based on rscaleFixed or rscaleCont or some mix of both?

[mattansb]

I've build a function to recreate the prior distributions of the BFmodel parameters - can be found here (BFEffects:::sample_from_priors).

If you get the chance to check this out, that would be amazing. Thanks!
(again, sorry for spamming)

[mattansb]

One more question regarding priors:
When modeling an interaction between a factor and a covariate (using generalTestBF), is the scale of the prior of the interaction term based on rscaleFixed or rscaleCont or some mix of both?

Seems that the rscale on interaction terms is a multiple of the rscale of the main effects:

library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************

data(puzzles)


set.seed(4)
result = anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID")


set.seed(4)
result2 = anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                  rscaleEffects = c("shape:color" = (sqrt(2)/2)^2))

set.seed(4)
result3 = anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                  rscaleEffects = c("shape:color" = sqrt(2)/2))

result[4]
#> Bayes factor analysis
#> --------------
#> [1] shape + color + shape:color + ID : 4.301216 ±2.34%
#> 
#> Against denominator:
#>   RT ~ ID 
#> ---
#> Bayes factor type: BFlinearModel, JZS
result2[4]
#> Bayes factor analysis
#> --------------
#> [1] shape + color + shape:color + ID : 4.301216 ±2.34%
#> 
#> Against denominator:
#>   RT ~ ID 
#> ---
#> Bayes factor type: BFlinearModel, JZS
result3[4]
#> Bayes factor analysis
#> --------------
#> [1] shape + color + shape:color + ID : 3.307638 ±2.37%
#> 
#> Against denominator:
#>   RT ~ ID 
#> ---
#> Bayes factor type: BFlinearModel, JZS

^{Created on 2019-05-23 by the reprex package (v0.2.1)}

[richarddmorey]

The prior on an interaction is the same as the prior type it becomes: e.g., fixed by fixed interaction is fixed, so it gets the fixed prior; a fixed by random interaction is random, so it gets the random prior; a fixed by continuous interaction is continuous, so it gets the continuous prior.

[mattansb]

Thanks for the reply!

If I understand correctly for what you say and from the following simulation:

1. An interaction of fixed by fixed does not gets a prior of fixed, but gets a prior of fixed^2

library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
#> 
#> Type BFManual() to open the manual.
#> ************

set.seed(4)
n <- 40
X <- rnorm(n)
X2 <- rnorm(n)
G <- rep(c(0,1), each = n/2)
G2 <- rep(c(0,1), times = n/2)
Y <- (X + G2)*(G + X2) + rnorm(n, sd = 0.5)

df <- data.frame(X = X,
                 X2 = X2,
                 G = factor(G),
                 G2 = factor(G2),
                 Y = Y)

# fixed by fixed ----------------------------------------------------------

set.seed(4)
BF1 <- lmBF(Y ~ G2*G, df, progress = FALSE)

set.seed(4)
BF2 <- lmBF(Y ~ G2*G, df, progress = FALSE,
            rscaleEffects = c("G2:G" = (sqrt(2)/2) * (sqrt(2)/2)))

set.seed(4)
BF3 <- lmBF(Y ~ G2*G, df, progress = FALSE,
            rscaleEffects = c("G2:G" = (sqrt(2)/2)))

extractBF(BF1)$bf # default
#> [1] 3.799659
extractBF(BF2)$bf # same as default
#> [1] 3.799659
extractBF(BF3)$bf
#> [1] 3.794392

2. An interaction of fixed by continuous gets a prior of a continuous

# fixed by cont -----------------------------------------------------------

set.seed(4)
BF1 <- lmBF(Y ~ X*G, df, progress = FALSE)

set.seed(4)
BF2 <- lmBF(Y ~ X*G, df, progress = FALSE,
            rscaleCont = c("X:G" = (sqrt(2)/4) * (sqrt(2)/2)))

set.seed(4)
BF3 <- lmBF(Y ~ X*G, df, progress = FALSE,
            rscaleCont = c("X:G" = (sqrt(2)/4)))

extractBF(BF1)$bf # default
#> [1] 6.18353
extractBF(BF2)$bf
#> [1] 7.229012
extractBF(BF3)$bf # same as default
#> [1] 6.18353

3. An interaction of continuous by continuous gets a prior of a continuous

# cont by cont ------------------------------------------------------------

set.seed(4)
BF1 <- lmBF(Y ~ X2*X, df, progress = FALSE)

set.seed(4)
BF2 <- lmBF(Y ~ X2*X, df, progress = FALSE,
            rscaleCont = c("X2:X" = (sqrt(2)/4) * (sqrt(2)/4)))

set.seed(4)
BF3 <- lmBF(Y ~ X2*X, df, progress = FALSE,
            rscaleCont = c("X2:X" = (sqrt(2)/4)))

extractBF(BF1)$bf # default
#> [1] 7850662
extractBF(BF2)$bf 
#> [1] 3141049
extractBF(BF3)$bf # same as default
#> [1] 7850662

^{Created on 2019-06-02 by the reprex package (v0.3.0)}

[richarddmorey]

The default scale for fixed factors is 1/2 (see

BayesFactor/pkg/BayesFactor/R/common.R

Lines 202 to 206 in 87013ba

	fixed = switch(priorType,
	ultrawide=1,
	wide=sqrt(2)/2,
	medium=1/2,
	stop("Unknown prior type.")),

) so it isn't fixed^2, it is just fixed.

[mattansb]

Oh! Thanks as always for your patience, Richard!