Regularisierte bayesianische logistische Regression in JAGS

13

Es gibt mehrere mathematisch anspruchsvolle Artikel, die das Bayes'sche Lasso beschreiben, aber ich möchte getesteten, korrekten JAGS-Code, den ich verwenden kann.

Könnte jemand einen Beispiel-BUGS / JAGS-Code veröffentlichen, der eine regulierte logistische Regression implementiert? Jedes Schema (L1, L2, Elasticnet) wäre toll, aber Lasso wird bevorzugt. Ich frage mich auch, ob es interessante alternative Implementierungsstrategien gibt.

Jack Tanner
quelle

Antworten:

19

Da die L1-Regularisierung einem Laplace (doppelt exponentiell) vor den relevanten Koeffizienten entspricht, können Sie dies wie folgt tun. Hier habe ich drei unabhängige Variablen x1, x2 und x3 und y ist die binäre Zielvariable. Auswahl des Regularisierungsparameters λ erfolgt hier durch Setzen eines Hyperpriors, in diesem Fall gerade gleichmäßig über einen großen Bereich.

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
}

Probieren wir es mit dem dclone Paket in R aus!

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)

Und hier sind die Ergebnisse im Vergleich zu einer unregelmäßigen logistischen 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 >>

Und wir können sehen, dass die drei bParameter tatsächlich gegen Null geschrumpft sind.

Ich weiß leider nicht viel über Prioritäten für den Hyperparameter der Laplace-Verteilung / den Regularisierungsparameter. Ich neige dazu, gleichmäßige Verteilungen zu verwenden und schaue auf den hinteren Teil, um zu sehen, ob er einigermaßen gut aussieht, z. Bisher war dies normalerweise der Fall. Es funktioniert auch für mich, es als Varianzparameter zu behandeln und die Empfehlungen von Gelman Prior-Verteilungen für Varianzparameter in hierarchischen Modellen zu verwenden.

Bogenschütze
quelle
1
Du bist der beste! Ich werde die Frage für eine Weile offen lassen, falls jemand eine andere Implementierung hat. Zum einen scheint es, dass binäre Indikatoren verwendet werden können, um variable Ein- / Ausschlüsse zu erzwingen. Dies kompensiert die Tatsache, dass unter Bayesian Lasso die Variablenauswahl nicht tatsächlich erfolgt, da die Betas mit dem doppelten Exponentialprior keine Posterioren haben, die genau Null sind.
Jack Tanner
Richtig, das mache ich auch. Es ist ähnlich wie beim Erstellen eines Priores mit einer Punktmasse von 0 und anschließender Verwendung des Nullentricks zum Abtasten des Prioresbichwird dann eine Mischung aus einer Punktmasse bei 0 und einem Dublettenexponential, obwohl der Code unterschiedlich ist. Ich bin gespannt, was andere Leute tun, +1 auf die Frage.
Bogenschütze