[R-sig-ME] Bayesian hypothesis testing using MCMCglmm
Alberto Gallano
alberto.gc8 at gmail.com
Mon Dec 22 01:27:17 CET 2014
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]]
More information about the R-sig-mixed-models
mailing list