Binomial Bayesian Regression with Metropolis Hastings in R

Question

I'm trying to implement the Metropolis Hastings algorithm in this problem but I'm having problems with the convergence.

$$Y_i|\beta_0,\beta_1 \sim \text{Binomial}(m_i,\theta_i)$$

where $logit(\theta) = \beta_0 + \beta_1x_i$.

$Y_i$ is the number of damage incidents of 6 possible (i.e number of successes). This also means $m_i = 6$.

I assume a normal distribution for the regression coefficients and I'm using random walk as a candidate distribution, so

$$\beta_0 \sim N(0,10)$$ $$\beta_0 \sim N(0,10)$$

$$\beta_j^*|\beta_{j}^{t-1} \sim N(\beta_{j}^{t-1},\sigma^2)$$

For the data I'm using

temp <- c(53,57,58,63,66,67,67,67,68,69,70,70,
          70,70,72,73,75,75,76,76,78,79,81)

success <- c(5,1,1,1,0,0,0,0,0,0,1,
             0,1,0,0,0,0,1,0,0,0,0,0)

failures <- c(1,5,5,5,6,6,6,6,6,6,
              5,6,5,6,6,6,6,5,6,6,6,6,6)

dat <- data.frame(success=success,failures=failures,temp=temp)

# Posterior distribution
post_beta <- function(Y,x,beta,m=6,mu_beta0=0,s2_beta0=10,mu_beta1=0,s2_beta1=10){
  b0 <- beta[1]
  b1 <- beta[2]
  pred <- b0+b1*x
  theta <- exp(pred)/(1+exp(pred))
  like <- dbinom(Y,m,prob = theta)
  prior <- dnorm(x=b0,mean = mu_beta0,sd = s2_beta0) *
    dnorm(x=b1,mean = mu_beta1,sd = s2_beta1)
  return(like*prior)
}

# Metropolis-Hastings

mh <- function(S,Y,x,init,fixed_sd){
  samples <- matrix(NA,S,2)
  colnames(samples) <- c("b0","b1")
  beta <- init
  for(s in 1:S){
    for(i in 1:length(beta)){
      can <- beta
      can[i] <- rnorm(1,mean = beta[i],fixed_sd)
      R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta))
      if(runif(1) <= R){
        beta <- can
      }
    }
    samples[s,] <- c(beta)
  }
  return(samples)
}

S <- 50000 # Simulations
Y <- dat$success
x <- dat$temp
init <- c(11.66299,-0.21623) # Initial values (values from glm output)
fixed_sd <- 0.01 # Fixed sd for the M-H

res <- mh(S = S,Y = Y,x = x,init = init,fixed_sd = fixed_sd)
plot(res[,1],type = "l",main="b0")
plot(res[,2],type = "l",main="b1")

The results using a frequentist framework using glm are the following.

res_freq = glm(cbind(data_p2$success,data_p2$failures)~x,family = binomial())
summary(res_freq )

And I tried to use this result as initial values but that doesn't seem to work to either. I have the feeling that the mistake could be in post_beta function.

Thanks in advance!

Answer 1

You can get nice grassy oscillation in the time series plots by changing the prior to

$$\begin{split}\beta_0&\sim \mathcal N(11.66299, .1^2)\\ \beta_1 &\sim \mathcal N(-0.21623, .1^2)\end{split}$$

and changing the standard deviation of your proposal/candidate to $.1$ for $\beta_0$ (keeping it at $.01$ for $\beta_1)$.

Lastly, the acceptance probability is missing the factor involving the candidate distribution.

temp <- c(53,57,58,63,66,67,67,67,68,69,70,70,
          70,70,72,73,75,75,76,76,78,79,81)

success <- c(5,1,1,1,0,0,0,0,0,0,1,
             0,1,0,0,0,0,1,0,0,0,0,0)

failures <- c(1,5,5,5,6,6,6,6,6,6,
              5,6,5,6,6,6,6,5,6,6,6,6,6)

dat <- data.frame(success=success,failures=failures,temp=temp)
# Posterior distribution
post_beta <- function(Y,x,beta,m=6,mu_beta0=11.66299,s2_beta0=.1,mu_beta1=-0.21623,s2_beta1=.1){
  b0 <- beta[1]
  b1 <- beta[2]
  pred <- b0+b1*x
  theta <- exp(pred)/(1+exp(pred))
  like <- dbinom(Y,m,prob = theta)
  prior <- dnorm(x=b0,mean = mu_beta0,sd = s2_beta0) *
    dnorm(x=b1,mean = mu_beta1,sd = s2_beta1)
  return(like*prior)
}

# Metropolis-Hastings

mh <- function(S,Y,x,init,fixed_sd){
  samples <- matrix(NA,S,2)
  colnames(samples) <- c("b0","b1")
  beta <- init
  for(s in 1:S){
    for(i in 1:length(beta)){
      can <- beta
      can[i] <- rnorm(1,mean = beta[i],fixed_sd[i])
      R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta)*dnorm(beta[i],can[i],fixed_sd[i])/dnorm(can[i],beta[i],fixed_sd[i]))
      if(runif(1) <= R){
        beta <- can
      }
    }
    samples[s,] <- c(beta)
  }
  return(samples)
}


S <- 50000 # Simulations
Y <- dat$success
x <- dat$temp
init <- c(11.66299,-0.21623) # Initial values (values from glm output)
fixed_sd <- c(.1, 0.01) # Fixed sd for the M-H

res <- mh(S = S,Y = Y,x = x,init = init,fixed_sd = fixed_sd)
par(mfrow=c(2,2))
plot(res[,1],type = "l",main="b0")
plot(res[,2],type = "l",main="b1")
acf(res[,1])
acf(res[,2])