Solved – Marginal effects from Bayesian probit

bayesianmarkov-chain-montecarloprobit

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.

Best Answer

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)
Related Question