Fitting a logistic regression

This example describes how to perform a Bayesian logistic regression using the malt function from the malt library. For an explanation on how to use the malt function see different vignette.

library(malt)

We are given \(n\) data-points which have \(d\) real-valued explanatory variables, summarised in a matrix \(X\in\mathbb{R}^{n\times d}\) and a vector of response variables \(Y\in\{0,1\}^{n}\). We want to find the \(\theta\in\mathbb{R}^d\) such that \((1+\exp(-\theta^\top X_i))^{-1}\) is a good predictor of the probability of \(Y_i=1\). We pick \(X\) and a true \(\theta\) and generate synthetic data \(Y\), which we then try to recover:

d=2
n=1000
X=matrix(rnorm(n*d),nr=n,nc=d)
true_theta=rep(1,d)
p=1/(1+exp(-X%*%true_theta))
Y=as.numeric(runif(length(p))<p)

The log-likelihood corresponding to the logistic regression is: \[ \log(L(y|\theta,x))~=~Y^\top X\theta~-~\sum_{i=1}^n \log\left(1+\exp(\theta^\top X_i)\right)\,. \] Together with a prior \(p(\theta)\) this generates a posterior potential \[ U(\theta)~=~-Y^\top X\theta~+~\sum_{i=1}^n \log\left(1+\exp(\theta^\top X_i)\right)~-~\log(p(\theta))\,. \] We are going to implement this in R with an impropper prior (we pre-compute \(Y^\top X\) to improve efficiency):

YX=c(Y%*%X)
U=function(theta){
-sum(YX*theta)+sum(log(1+exp(X%*%theta)))
}
grad=function(theta){
-YX+c(c(1/(1+exp(-X%*%theta)))%*%X)
}
init=rep(0,d)
n_steps=1000
g=1
h=0.1
L=floor(2/h)
output=malt(init, U, grad, n_steps, g, h, L)
chain=output$samples
apply(chain,2,mean)
#> [1] 0.9575116 1.0077304
output$acceptance
#> NULL
plot(chain[,1],type="l")

plot(chain[,2],type="l")