Priors on odds and probability of success

In Bayesian Approaches to Clinical Trials and Health-Care Evaluation (David J. Spiegelhalter, Keith R. Abrams, Jonathan P. Myles) they mention that an non informative prior should be uniform over the range of interest. They combine this with the desire that Odds ratio has the same prior as 1/Odds ratio and from this they conclude log Odds ratio should be uniform distributed. They write it is equivalent to Beta(0,0), hence not a proper distribution, but use it anyway. In this post looks at the prior both at log odds scale and probability scale.
As a side, this exercise forced me to look at the Jacobian, I needed it to transform a density to the other scale. For this yacas came in handy, it has been a while since I differentiated by hand.

Analysis of prior

In general, I will try to give both the distribution according to a MCMC sampler (black lines) and an analytical solution (green line). The MCMC sampler used is LaplacesDemon, since it is fully in R it allows reusing analytical code and easy adding of the Jacobian.
As a side, in many of the plots I made the smoothing bandwidth a bit smaller, so some of the details could be captured.
Four priors are used. Probability of success Beta(0.5,0.5) and Beta(1,1). Log odds uniform on -10 to 10 and normal(0,10).
I think LaplacesDemon with the standard algorithm has some difficulty with the beta(0.5,0.5) distribution.

Using Data

After observing data, one success and five failure again same prior. For the priors on log odds there is no exact formulas to obtain the exact likelihood in the book. I used a formula for normal approximation of log odds ratio (2.30) and simplified from there. I am actually surprised how close that is to the MCMC sampler.

Code

# general code
library(magrittr)
library('LaplacesDemon')

#conversion and Jacobian
prob2Lo <- function(prob) log(prob/(1-prob))
Lo2prob <- function(Lo) exp(Lo)/(1+exp(Lo))
Jprob2Lo <- function(prob) prob*(1-prob)
JLo2prob <- function(Lo) (exp(Lo)+1)^2/exp(Lo)

# wrapper for Jacobian
dprob2dLo <- function(Lo,dprob,log=FALSE) {
  prob <- Lo2prob(Lo)
  if (log) {dprob(prob)+log(Jprob2Lo(prob))
  } else dprob(prob)*Jprob2Lo(prob)
}

dLo2dprob <- function(prob,dLo,log=FALSE) {
  Lo <- prob2Lo(prob)
  if (log) {dLo(Lo)+log(JLo2prob(Lo))
  } else dLo(Lo)*JLo2prob(Lo)
}

# helper function
densplot <- function (x,adjust=1,...) {
  density(x,adjust) %>%
      plot(.,...)
}

mon.names <- "LP"
parm.names <- c('prob')
MyData <- list(mon.names=mon.names,
    parm.names=parm.names)
N <-1
Model <- function(parm, Data)
{
  LL <- 1
  yhat <- 1
  LP <- dbeta(parm,.5,.5,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(0.5)
Fit <- LaplacesDemon(Model, Data=MyData,
    Initial.Values,
    Iteration=1e5,Status=1e4)
png('B55.png')
par(mfrow=c(1,2))
densplot(prob2Lo(Fit$Posterior2),
    adjust=.2,
    main='Prior: Probability ~ Beta .5 .5',
    xlab='Log odds')
Lo <- seq(-10,10,.1)
dLo <- dprob2dLo(Lo,dprob=function(prob)
      dbeta(prob,1,1))
lines(x=Lo,y=dLo,col='green')
##
densplot(Fit$Posterior2,
    adjust=.01,
    main='Prior: Probability ~ Beta .5 .5',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dbeta(prob,.5,.5)
lines(x=prob,y=dprob,col='green')
dev.off()
###
Model <- function(parm, Data)
{
  LL <- 1
  yhat <- 1
  LP <- dbeta(parm,1,1,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(0.5)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values,
    Iteration=1e5,Status=1e4)
png('B11.png')
par(mfrow=c(1,2))
densplot(prob2Lo(Fit$Posterior2),
    adjust=.1,
    main='Prior: Prob ~ Beta 1 1',
    xlab='Log odds')
Lo <- seq(-10,10,.1)
dLo <- dprob2dLo(Lo,dprob=function(prob)
      dbeta(prob,1,1))
lines(x=Lo,y=dLo,col='green')
densplot(Fit$Posterior2,adj=.01,
    main='Prior: Prob ~ Beta 1 1',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dbeta(prob,1,1)
lines(x=prob,y=dprob,col='green')
dev.off()
###
# prior on log odds
##
Model <- function(parm, Data)
{
  LL <- 1
  yhat <- 1
  LP <- dunif(parm,-10,10,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(0)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values,
    Iteration=1e5,
    Status=1e4)
png('u1010.png')
par(mfrow=c(1,2))
densplot(Fit$Posterior2,
    adjust=.1,
    main='Prior: Log Odds ~ Unif -10 10',
    xlab='Log odds')
Lo <- seq(-10.1,10.1,.1)
dLo <- dunif(Lo,-10,10)
lines(x=Lo,y=dLo,col='green')
densplot(Lo2prob(Fit$Posterior2),
    adjust=.01,
    main='Prior: Log Odds ~ Unif -10 10',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dLo2dprob(prob,dLo=function(Lo)
      dunif(Lo,-10,10))
lines(x=prob,y=dprob,col='green')
dev.off()
#
Model <- function(parm, Data)
{
  LL <- 1
  yhat <- 1
  LP <- dnorm(parm,0,10,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(0)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values,
    Iteration=1e5,
    Status=1e4)
png('N010.png')
par(mfrow=c(1,2))
densplot(Fit$Posterior2,
    adjust=.1,
    main='Prior: Log Odds ~ norm 0,10',
    xlab='Log odds')
Lo <- seq(-20.1,20.1,.1)
dLo <- dnorm(Lo,0,10)
lines(x=Lo,y=dLo,col='green')
##
densplot(Lo2prob(Fit$Posterior2),
    adjust=.01,
    main='Prior: Log Odds ~ norm 0,10',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dLo2dprob(prob,dLo=function(Lo) dnorm(Lo,0,10))
lines(x=prob,y=dprob,col='green')
dev.off()

#############
# with data
#############
# 1 out of 6
parm.names <- c('prob')
MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    n=6,s=1)

Model <- function(parm, Data)
{
  prob <- parm
  LL <- dbinom(prob,x=Data$s,size=Data$n,log=TRUE)
  LP <- dbeta(prob,.5,.5,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(.5)
Fit <- LaplacesDemon(Model, Data=MyData,
    Initial.Values,
    Iteration=1e5,Status=1e4)
png('db55.png')
par(mfrow=c(1,2))
densplot(prob2Lo(Fit$Posterior2),
    adjust=.1,
    main='Prior: Probability ~ Beta .5 .5',
    xlab='Log odds')
Lo <- seq(-10,10,.1)
dLo <- dprob2dLo(Lo,dprob=function(prob)
      dbeta(prob,1.5,5.5))
lines(x=Lo,y=dLo,col='green')
##
densplot(Fit$Posterior2,
    adjust=.01,
    main='Prior Probability ~ Beta .5 .5',
    xlab='Probability Success',
    xlim=c(0,1))
prob <- seq(0.01,.99,.01)
dprob <- dbeta(prob,1.5,5.5)
lines(x=prob,y=dprob,col='green')
dev.off()
# 1 1
parm.names <- c('prob')
MyData <- list(mon.names=mon.names,
    parm.names=parm.names,n=6,s=1)

Model <- function(parm, Data)
{
  prob <- parm
  LL <- dbinom(prob,x=Data$s,size=Data$n,log=TRUE)
  LP <- dbeta(prob,1,1,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(.5)
Fit <- LaplacesDemon(Model, Data=MyData,
    Initial.Values,
    Iteration=1e5,Status=1e4)
png('db11.png')
par(mfrow=c(1,2))
densplot(prob2Lo(Fit$Posterior2),
    adjust=.1,
    main='Prior: Probability ~ Beta 1 1',
    xlab='Log odds')
Lo <- seq(-10,10,.1)
dLo <- dprob2dLo(Lo,dprob=function(prob) dbeta(prob,2,6))
lines(x=Lo,y=dLo,col='green')
densplot(Fit$Posterior2,
    adjust=.01,
    main='Prior Probability ~ Beta 1 1',
    xlab='Probability Success',
    xlim=c(0,1))
prob <- seq(0.01,.99,.01)
dprob <- dbeta(prob,2,6)
lines(x=prob,y=dprob,col='green')
dev.off()
##u -10 10 #########
mon.names <- "LP"
parm.names <- c('Lo')
MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    n=6,s=1)

Model <- function(parm, Data)
{
  Lo <- parm
  LL <-dprob2dLo(Lo,
      dprob=function(prob)
        dbinom(prob,
            x=Data$s,
            size=Data$n ,
            log=TRUE),
      log=TRUE)
  yhat <- Lo
  LP <- dunif(Lo,-10,10,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(.1)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values,
    Iteration=1e5,
    Status=1e4)
png('du1010.png')
par(mfrow=c(1,2))
densplot(Fit$Posterior2,
    adjust=.1,
    main='Prior: Log Odds ~ Unif -10 10',
    xlab='Log odds')
Lo <- seq(-20.1,20.1,.1)
dLo <- dnorm(Lo,log(1.5/5.5),sqrt(1/1.5+1/5.5))
lines(x=Lo,y=dLo,col='green')
densplot(Lo2prob(Fit$Posterior2),
    adjust=.01,
    main='Prior: Log Odds ~ Unif -10 10',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dLo2dprob(prob,dLo=function(Lo) dnorm(Lo,log(1.5/5.5),sqrt(1/1.5+1/5.5)))
lines(x=prob,y=dprob,col='green')
dev.off()
## N 0 10
parm.names <- c('Lo')
MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    n=6,s=1)

Model <- function(parm, Data)
{
  Lo <- parm
  LL <-dprob2dLo(Lo,
      dprob=function(prob)
        dbinom(prob,
            x=Data$s,
            size=Data$n,
            log=TRUE),
      log=TRUE)
  yhat <- Lo
  LP <- dnorm(Lo,0,10,log=TRUE)+LL
  Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LL,
      yhat=1, parm=parm)
  return(Modelout)
}
Initial.Values <- c(.1)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values,
    Iteration=1e5,
    Status=1e4)
png('dn010.png')
par(mfrow=c(1,2))

t1 <- 0
t2 <- log(1.5/5.5)
var1 <- 100
var2 <- 1/1.5+1/5.5
tt <- (t1/var1+t2/var2)/(1/var1+1/var2)
vv <- 1/(1/var1+1/var2)
densplot(Fit$Posterior2,
    adjust=.1,
    main='Prior: Log Odds ~ norm 0,10',
    xlab='Log odds')
Lo <- seq(-20.1,20.1,.1)
dLo <- dnorm(Lo,tt,sqrt(vv))
lines(x=Lo,y=dLo,col='green')
densplot(Lo2prob(Fit$Posterior2),
    adjust=.01,
    main='Prior: Log Odds ~ norm 0,10',
    xlab='Probability Success')
prob <- seq(0.01,.99,.01)
dprob <- dLo2dprob(prob,dLo=function(Lo) dnorm(Lo,tt,sqrt(vv)))
lines(x=prob,y=dprob,col='green')
dev.off()