############################################################################
## Code and examples by Vanja Dukic, U of Colorado at Boulder
## Bayesian Statistics and Computing 2019
##
## LECTURE 7 - Priors in logistic regression
## Differences between the log odds and probability scale
############################################################################


# Logistic regression: logit(theta) = alpha0 + alpha1 X_i
# Let X_i be a binary indicator: treatment/placebo
# So two regressions: 
# Placebo group: logit(theta) = alpha0 
# Treatment group: logit(theta) = alpha0 + alpha1


# It is tempting to set priors on alpha0 and aplha1 as "vague"
# For example, that can be Normal with mean 0 and variance 100
# But alphas are on the log-odds scale, which is difficult to work with
# Easiest to think on a probability scale

# So think about what a prior on alpha0 implies for the prior on theta

# logit(theta) = alpha0   => theta = inverse-logit(alpha0):
# recall that logit(p) = log-odds, or log(p/(1-p))
# so if log(p/(1-p)) = alpha0 then p/(1-p) = exp(alpha0), and p = exp(alpha0)/(1+exp(alpha0))

# You can do a formal change of measure here, using the change of measure formula

# Below are some simulations to help you understand the impact of 
# different alpha priors on priors on theta 
# 


x11()
ralpha0 = rnorm(1000,0,100)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,100)')
hist(rtheta)



x11()
ralpha0 = rnorm(1000,0,10)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,10)')
hist(rtheta)



x11()
ralpha0 = rnorm(1000,0,4)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,4)')
hist(rtheta)



x11()
ralpha0 = rnorm(1000,0,2)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,2)')
hist(rtheta)


x11()
ralpha0 = rnorm(1000,0,1.5)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,1.5)')
hist(rtheta)



x11()
ralpha0 = rnorm(1000,0,1)
rtheta = exp(ralpha0)/(1+exp(ralpha0))
par(mfrow=c(1,2))
hist(ralpha0, main='N(0,1)')
hist(rtheta)