15
$\begingroup$

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.

$\endgroup$
20
$\begingroup$

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.adapt=1000, 
                 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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Jack Tanner May 17 '12 at 2:46
  • $\begingroup$ 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. $\endgroup$ – jbowman May 17 '12 at 3:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.