Marginal effects from Bayesian probit I'm trying to run a standard Bayesian probit model, and I can't find any packages in R that will give me marginal effects (the most common way to interpret probit results in my field), nor do they give me the elements I need to estimate them. According to Koop et al , page 209:

The posterior distribution of the marginal effect can be obtained by
  calculating and collecting the values
  $(\beta_j^{(i)}\phi(x\beta^{(i)}))^{1500}_{i=1}$, where $\beta^{(i)}$
  represents the ith postconvergence draw from the Gibbs sampler.

My interpretation of this is, if I can get the software to output the linear prediction (i.e. the latent variable in the latent variable interpretation of the probit model), for each draw, and multiple that by the coefficient on the variable of interest, I will get that variable's marginal effect on each observation in the draw. I can then aggregate within the draw to get the average effect for that draw, and aggregate those across all draws to get the posterior distribution of the average marginal effect.
I have so far not found a Bayesian package in R that can get me this. I could probably program WinBUGs to do so, but I'm trying to stick to an entirely R solution if possible (I work with a remote team and they know R but getting up on BUGs will be a big lift). With all the Bayesian packages out there for R these days, I feel like someone should have had a solution to this, and I just haven't discovered it.
 A: From this post, all you need to calculate the quantity above is the samples for $\beta$ and a value for x, the covariate vector. 
One option is the bayes.probit function from the LearnBayes package which provides posterior draws of the coefficients.
Here is an example, using the example in the help file for bayes.probit:
library(LearnBayes)
response=c(0,1,0,0,0,1,1,1,1,1)
covariate=c(1,2,3,4,5,6,7,8,9,10)
X=cbind(1,covariate)
prior=list(beta=c(0,0),P=diag(c(.5,10)))
m=1500
s=bayes.probit(response,X,m,prior)

x = c(1,1)
marginal_effect = dnorm(s$beta %*% x) * s$beta[,2]
mean(marginal_effect)

A: spatialprobit package in R does this with the function impacts(model) and is also nice because it implements sparse matrices for calculation.  I have found it straightforward to use.  The impacts function will calculate direct, indirect, and total impacts.  This is all if you are doing a spatial probit AR model, but I would imagine you could set it up to work without any spatial neighbors quite easily.  Anyway, thought I would add this for people looking for spatial bayes answers to the same problem.
