# Regularized bayesian logistic regression in JAGS

There are several math-heavy papers that describe the Bayesian Lasso, but I want tested, correct JAGS code that I can use.

Could someone post sample BUGS / JAGS code that implements regularized logistic regression? Any scheme (L1, L2, Elasticnet) would be great, but Lasso is preferred. I also wonder if there are interesting alternative implementation strategies.

Since L1 regularization is equivalent to a Laplace (double exponential) prior on the relevant coefficients, you can do it as follows. Here I have three independent variables x1, x2, and x3, and y is the binary target variable. Selection of the regularization parameter $\lambda$ is done here by putting a hyperprior on it, in this case just uniform over a good-sized range.

model {
# Likelihood
for (i in 1:N) {
y[i] ~ dbern(p[i])

logit(p[i]) <- b0 + b[1]*x1[i] + b[2]*x2[i] + b[3]*x3[i]
}

# Prior on constant term
b0 ~ dnorm(0,0.1)

# L1 regularization == a Laplace (double exponential) prior
for (j in 1:3) {
b[j] ~ ddexp(0, lambda)
}

lambda ~ dunif(0.001,10)
# Alternatively, specify lambda via lambda <- 1 or some such
}


Let's try it out using the dclone package in R!

library(dclone)

x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)

prob <- exp(x1+x2+x3) / (1+exp(x1+x2+x3))
y <- rbinom(100, 1, prob)

data.list <- list(
y = y,
x1 = x1, x2 = x2, x3 = x3,
N = length(y)
)

params = c("b0", "b", "lambda")

temp <- jags.fit(data.list,
params=params,
model="modela.jags",
n.chains=3,
n.update=1000,
thin=10,
n.iter=10000)


And here are the results, compared to an unregularized logistic regression:

> summary(temp)

<< blah, blah, blah >>

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean     SD Naive SE Time-series SE
b[1]   1.21064 0.3279 0.005987       0.005641
b[2]   0.64730 0.3192 0.005827       0.006014
b[3]   1.25340 0.3217 0.005873       0.006357
b0     0.03313 0.2497 0.004558       0.005580
lambda 1.34334 0.7851 0.014333       0.014999

2. Quantiles for each variable: << deleted to save space >>

> summary(glm(y~x1+x2+x3, family="binomial"))

<< blah, blah, blah >>

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.02784    0.25832   0.108   0.9142
x1           1.34955    0.32845   4.109 3.98e-05 ***
x2           0.78031    0.32191   2.424   0.0154 *
x3           1.39065    0.32863   4.232 2.32e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

<< more stuff deleted to save space >>


And we can see that the three b parameters have indeed been shrunk towards zero.

I don't know much about priors for the hyperparameter of the Laplace distribution / the regularization parameter, I'm sorry to say. I tend to use uniform distributions and look at the posterior to see if it looks reasonably well-behaved, e.g., not piled up near an endpoint and pretty much peaked in the middle w/o horrible skewness problems. So far, that's typically been the case. Treating it as a variance parameter and using the recommendation(s) by Gelman Prior distributions for variance parameters in hierarchical models works for me, too.

• You're the best! I'll leave the question open for a little while in case somebody has another implementation. For one, it seems that binary indicators can be used to impose variable inclusion / exclusion. This compensates for the fact that under Bayesian Lasso variable selection doesn't actually happen, since the betas with the double exponential prior will not have posteriors that are exactly zero. – Jack Tanner May 17 '12 at 2:46
• Right, I do that too. It's similar to creating a prior with a point mass at 0, then using the zeroes trick to sample from it (the prior on $b_i$ then becomes a mixture of a point mass at 0 and a duble exponential), although the code is different. I'll be interested to see what other people do, +1 to the question. – jbowman May 17 '12 at 3:01