Uregulowana bayesowska regresja logistyczna w JAGS

Istnieje kilka prac matematycznych opisujących lasso bayesowskie, ale chcę przetestować poprawny kod JAGS, którego mogę użyć.

Czy ktoś może opublikować próbkę kodu BŁĘDY / JAGS, który implementuje regulowaną regresję logistyczną? Każdy schemat (L1, L2, Elasticnet) byłby świetny, ale preferowany jest Lasso. Zastanawiam się także, czy istnieją ciekawe alternatywne strategie wdrażania.

Ponieważ regularyzacja L1 jest równoważna Laplace'owi (podwójnemu wykładniczemu) przed odpowiednimi współczynnikami, możesz to zrobić w następujący sposób. Tutaj mam trzy niezależne zmienne x1, x2 i x3, a y jest binarną zmienną docelową. Wyboru parametru regularyzacji dokonuje się tutaj poprzez nałożenie na niego hiperpriora, w tym przypadku po prostu jednolity w sporej wielkości.$\lambda$

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
}

Wypróbujmy to, używając dclonepakietu w 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)

A oto wyniki, w porównaniu do nieregularnej regresji logistycznej:

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

I widzimy, że trzy bparametry rzeczywiście zostały zmniejszone do zera.

Przykro mi to mówić, ale niewiele wiem o priorytetach hiperparametru rozkładu Laplace'a / parametru regularyzacji. Zwykle używam równomiernych rozkładów i patrzę na tył, aby zobaczyć, czy wygląda dość dobrze, np. Nie jest ułożony w pobliżu punktu końcowego i prawie osiąga szczyt w środkowej części bez okropnych problemów ze skośnością. Jak dotąd zwykle tak jest. Traktowanie go jako parametru wariancji i korzystanie z rekomendacji Gelmana Wcześniejsze rozkłady parametrów wariancji w modelach hierarchicznych również mi się podobają.