[R-sig-ME] Bayesian hypothesis testing using MCMCglmm

[Alberto Gallano]

I'm fitting a linear mixed model using MCMCglmm and I want to determine the
probability that some a prior hypotheses I have about the slope and
intercept of the fixed effects are true. The slope estimate from my fitted
model is 1.21 (HDI 1.04, 1.38), while the intercept is -0.68 (HDI -0.46,
-0.93). I want to calculate the probability of getting a slope of 1 and an
intercept of -0.5.

One of the often touted main differences between frequentist and Bayesian
frameworks is the latter's ability to quantify the probability of a
hypothesis, given data. From what i've read, there seems to be two main
approaches to operationalizing this:

1) model comparison, where one model is fitted with a very narrow,
"spiked", prior for the hypothesis under consideration and another model
fitted with a much broader prior. The DICs from each model are presumably
then compared to determine which is better. I don't really understand where
the probability for the hypothesis comes from here.

2) determine a Region of Practical Equivalence (ROPE) around the hypothesis
to be tested, then see if the 95% HDI is encompassed within this region or
vice versa. Again, I don't really understand where the probability for the
hypothesis comes from here.

If I were doing this in a frequentist framework, i'd see if the 95% CI of
my estimates encompassed the hypothesized slope/intercept, and use that as
a hypothesis test. In a Bayesian framework, I should (I think) be able to
do better than that and actually get a probability for the hypothesized
slope/intercept, given the data.

My question is, how would I do this in practice, using MCMCglmm?

My model is below:

priors1 <- list(
    B = list(mu = rep(0, 3), V = diag(9, 3)),
    G = list(G1 = list(V = 1, nu = 0.002), G2 = list(V = diag(2)/2, nu =
0.002)),
    R = list(V = 1, nu = 0.002)
    )

# inverse of sigma matrix of phylogenetic correlation
inv_phylo_mat <- inverseA(tree, nodes = "TIPS", scale = TRUE)

fit <- MCMCglmm(
    fixed = I1.M1 ~ I2.M1.species.mean + I2.M1.within.species,
    rcov = ~ units,
    random = ~ phylo + idh(1+I2.M1.within.species):species,
    data = incisor.dat,
    family = "gaussian",
    ginverse = list(phylo = inv_phylo_mat$Ainv),
    prior = priors1,
    pr = TRUE,
    pl = TRUE,
    nitt = 1.1e+6, thin = 10, burnin = 1e+5,
    verbose = FALSE
    )

best,
Alberto

	[[alternative HTML version deleted]]