# How to set the gamma prior around the MLE estimate?

I have a MLE estimate for the lambda of the possion distribution.

However, the size of my data sets is very small. I choose use the gamma prior in jags. I want to set my prior around the MLE estimate. For example MLE estimate is 5. How can I set my prior using dgamma? Especially for the variance I understand that this is the method called Empirical Bayes.

The mean and variance of a gamma distribution (beware of the parametrizations, see https://en.wikipedia.org/wiki/Gamma_distribution) are $$k\theta$$ and $$k\theta^2$$. Hence, if you set $$\theta=5/k$$ you will get a prior mean of five. The higher you set $$k$$, the more the prior will be centered around 5, as the variance becomes $$k\cdot\frac{25}{k^2}=\frac{25}{k}$$ Alternatively, "setting my prior around the MLE estimate" could be taken to imply to set it around the mode, which is given by $$(k-1)\theta$$. Hence, setting $$\theta=5/(k-1)$$ would be the way to go here, leading to a variance of $$k\cdot\frac{25}{(k-1)^2},$$ which however qualitatively leads to the same conclusions.

Yet, I would not recommend this procedure. Empirical Bayes usually proceeds on borrowing evidence from related experiments rather than tweak the prior to make the posterior look very centered around the MLE.

Here, your readers might be tempted to conclude that you wanted to see the result 5 in the first place. If you, a priori, did not have strong reasons to believe in a point estimate 5 (if you had a very strong reason, you might skip analysing the data altogether...), why do you suddenly have such reasons after looking at the MLE that, as you point out, results from a sample that, due to its size, you cannnot trust too much?

You can play around with the following code to see the effect of modifying k:

x <- seq(0,10,by=0.1)
k <- c(30, 15, 3)
theta <- 5/k
cols <- c("darkgreen", "lightblue", "salmon")
plot(x, dgamma(x, shape = k[1], scale = theta[1]), type="l", lwd=2, col=cols[1])
lines(x, dgamma(x, shape = k[2], scale = theta[2]), type="l", lwd=2, col=cols[2])
lines(x, dgamma(x, shape = k[3], scale = theta[3]), type="l", lwd=2, col=cols[3])
abline(v=5, lty=2, lwd=2, col="brown")
abline(v=(k-1)*theta, lty=2, col=cols, lwd=2) # the modes
legend("topright", legend=c(paste0("k=",k), "mean"), col=c(cols, "brown"), lty=2, lwd=2)