I have a philosophical question about interpreting samples from a posterior. Consider the simple example:
library(rstanarm)
set.seed(1839)
x <- rnorm(100)
y <- x + rnorm(100, 0, 4)
dat <- data.frame(x, y)
model <- stan_glm(
y ~ x,
prior = normal(0, 2),
prior_intercept = normal(0, 2),
prior_aux = cauchy(0, 5)
)
Let’s say my client wants to know if x
positively predicts y
. I could calculate how many of my posterior samples were above zero:
mean(as.matrix(model)[, "x"] > 0)
Which returns 0.997
.
Which would you all consider to be an accurate characterization of the probability? Consider the interpretations:
A) “The probability that there is a positive relationship between x and y is .997”
B) “The probability that there is a positive relationship between x and y, given the data, is .997”
C) “The posterior probability that there is a positive relationship between x and y, is .997” (Problem here is that clients won’t know what this means, likely)
D) “The probability that there is a positive relationship between x and y, given our model and priors, is .997”
E) “Given the model and priors we specified, as well as the data we observed, there is a .997 probability that there is a positive relationship between x and y”
(Nevermind for now that an estimate of .0001 might as well be zero, but is still considered “positive” for the time being).
I really enjoy MCMC and sampling from the posterior, because it leads to intuitive probability statements. However, I do not want to overstate the results, either. What would you consider to be the proper way to interpret proportions form a posterior in probabilistic terms?
Other ways you could look at the data, which I assume would follow the same probability statements:
> mean(as.matrix(model)[, "x"] > 0.5)
[1] 0.9415
> mean(as.matrix(model)[, "x"] > 0.5 & as.matrix(model)[, "x"] < 1.5)
[1] 0.7185
> mean(as.matrix(model)[, "x"] < 0.5 & as.matrix(model)[, "x"] < 1.5)
[1] 0.0585