The arm
package includes the sim()
function. On its R help page it states:
library(arm)
?sim
"This generic function gets posterior simulations of sigma and beta
from a ‘lm’ object, or simulations of beta from a ‘glm’ object, or
simulations of beta from a ‘merMod’ object"
From what I gathered sim()
does not necessarily use any form of MCMC. It only uses MCMC when the posterior distribution cannot be defined analytically. Now, given some basic knowledge about MCMC: It uses the values for the initial parameters, the given data, the model and the prior distribution to determine whether the initial values for the parameters are good values or not. The MCMC algorithm will (ideally) converge towards the joint posterior distribution of the parameters given the priors, the initial parameter values, the model and the data. Thereby fitting the model, and sampling from the posterior distribution at the same time.
The fact that sim()
only uses MCMC when the posterior distribution cannot be defined analytically seems to imply that sim()
uses maximum a posteriori probability estimation (MAP) in all cases. This, in my opinion, has the following reason. I assume that sim()
uses non-informative priors in all cases. If sim()
uses MAP this has the effect that the problem comes down to maximizing the likelihood function. Which in turn entails that there will be a marginal numerical difference between e.g. confidence intervals and credible intervals. Computationally this would be a lot more tractable than EAP (expected a posteriori estimation).
But how does sim()
come to a suitable interval that serves as a non-informative prior for a given model? Is there a simple mathematical method that allows one to do this for any given model?
Don’t hesitate to include complicated math if it need be!
(This is probably a post that could be on both CrossValidated and Stackoverflow. If there are strong opinions about this I will flag my own post and ask for it to be moved. Just tell me in the comments.)