Régression logistique bayésienne régularisée dans JAGS

Puisque la régularisation L1 est équivalente à un a priori de Laplace (double exponentielle) sur les coefficients pertinents, vous pouvez le faire comme suit. Ici, j'ai trois variables indépendantes x1, x2 et x3, et y est la variable cible binaire. La sélection du paramètre de régularisation se fait ici en y mettant un hyperprior, dans ce cas juste uniforme sur une plage de bonne taille. $\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
}

Essayons-le en utilisant le dclonepackage dans 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)

Et voici les résultats, par rapport à une régression logistique non régularisée:

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

Et nous pouvons voir que les trois bparamètres ont en effet été réduits à zéro.

Je ne connais pas grand chose aux priors pour l'hyperparamètre de la distribution de Laplace / le paramètre de régularisation, je suis désolé de le dire. J'ai tendance à utiliser des distributions uniformes et à regarder la partie postérieure pour voir si elle semble raisonnablement bien se comporter, par exemple, non empilée près d'un point final et à peu près au maximum au milieu sans horribles problèmes d'asymétrie. Jusqu'à présent, cela a généralement été le cas. Le traiter comme un paramètre de variance et utiliser la ou les recommandations de Gelman Prior pour les paramètres de variance dans les modèles hiérarchiques fonctionne aussi pour moi.

jbowman
la source

Tu es le meilleur! Je vais laisser la question ouverte pendant un petit moment au cas où quelqu'un aurait une autre implémentation. D'une part, il semble que les indicateurs binaires peuvent être utilisés pour imposer une inclusion / exclusion variable. Cela compense le fait que, sous le lasso bayésien, la sélection des variables ne se produit pas réellement, car les bêtas avec le double exponentiel avant n'auront pas de postérieurs qui sont exactement zéro.

Jack Tanner

b_{i}

$b_i$ devient alors un mélange d'une masse ponctuelle à 0 et d'une exponentielle double), bien que le code soit différent. Je serai intéressé de voir ce que font les autres, +1 à la question.

jbowman

Régression logistique bayésienne régularisée dans JAGS

Réponses: