R implementation of Polson 2013 "Bayesian inference for logistic models using Polya-Gamma latent variables"

gibbs_logistic.R

library(BayesLogit)
library(ggplot2)
#y: response variable
#n: binomial size parameter
#X: explanatory design matrix
#lambda: prior parameter
gibbs_logistic <- function(y,n,X,lambda,iter){
  N <- length(y)
  d <- ncol(X)
  B <- diag(lambda, d)
  beta_tilde <- matrix(0, iter, d)
  beta_tilde[1,] <- rnorm(d)
  kappa <- (y-0.5*n)
  for(i in 2:iter){
    eta <- X%*%beta_tilde[i-1,]
    Omega <- diag(rpg(N, n, eta))
    Vinv <- (t(X)%*%Omega%*%X + B)
    U <- chol(Vinv)
    A <- forwardsolve(t(U), t(X)%*%(kappa))
    m <- backsolve(U, A) #equivalent to #m <- solve(Vinv%*%(t(X)%*%(y-0.5*n)))
    beta_tilde[i,]  = m + backsolve(U, rnorm(d)) #multiply to inverse of U
  }
  return(beta_tilde)
}

set.seed(123)
beta <- c(1,1)
X <- matrix(runif(2*100),100,2)
y <- rbinom(100,10,plogis(X%*%beta))

beta_sample <- gibbs_logistic(y,10,X,1,2000)

dfs <- expand.grid(row=1:2,iter=1:2000)
dfs$value <- as.vector(beta_sample)

ggplot(dfs, aes(x=iter, y=value))+
  geom_line(colour="grey")+
  facet_grid(row~., scales = "free", labeller = label_both)+
  theme_classic(14)+
  theme(strip.text.y = element_text(angle = 0),
        axis.text = element_text(colour = "black"))

lp <- function(beta){
  sum(dbinom(y,10,plogis(X%*%beta), log = TRUE))+sum(dnorm(beta,log=TRUE))  
}

df <- expand.grid(b1=seq(0.4,1.5,by=0.02),b2=seq(0.4,1.6,by=0.02))
lp_v <- apply(df, 1, function(b){lp(c(b[1],b[2]))})

burnin <- 1:500
dfs <- data.frame(beta_sample[-burnin,])
colnames(dfs) <- c("b1","b2")
ggplot(dfs, aes(x=b1, y=b2))+
  geom_point(alpha=0.2)+
  geom_contour(data = df, aes(z=exp(lp_v), colour=after_stat(level)))+
  scale_color_viridis_c()+
  theme_classic(16)

#ggsave("post.png")