# Binomial Bayesian Regression with Metropolis Hastings in R

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.

• this line R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta)) is missing a correction factor I think Jun 24 at 2:40
• post betta doesn't need the prior involved either because it cancels out in the acceptance probability. that being said, convergence is tough for these kinds of problems Jun 24 at 3:26
• @Jellyfish thanks for the reply! What I don't know is if the likelihood is well defined in my post_beta function and the other thing is perhaps the standard deviation I'm using for the candidate is not the right one.
– Seb
Jun 24 at 5:45

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])


• Thanks! On the factor involving the candidate distribution, I didn't consider it since I'm using a symmetric distribution for the candidate this line cancels out dnorm(beta[i],can[i],fixed_sd[i])/dnorm(can[i],beta[i],fixed_sd[i]). One more thing, how did you come up with those values for the standard deviation?
– Seb
Jun 24 at 17:05
• @Seb By trial and error. I chose priors close to the frequentist framework with small standard deviation. Then based on the magnitudes of the mean, I chose the standard deviation of the proposal. I think convergence might be easier if the model fit the data better, but I’m not sure, because using Bayesian software like JAGS would probably not require such a specific prior and still have good results. So, I’m not entirely sure either Jun 24 at 18:30