Bayesianische Variablenselektion

Variablenselektion

Oftmals hat man sehr viele Kovariablen und will nur die “relevanten” auswählen. Als Beispiel arbeiten wir hier mit folgenden Datensatz (Quelle):

Für 442 Diabetes-Patienten wurden ein Maß \(y_i\) für den Fortschritt der Krankheit sowie folgende Kovariablen bestimmt:

Alter (age)
Geschlecht (sex)
BMI (bmi)
Mittlerer Blutdruck (map)
Sechs verschiedene Blutwerte (tc, ldl, hdl, tch, ltg, glu)

Alle Variablen wurden standardisiert (Mittelwert 0, Standardabweichung 1)

summary(y)

##        V1         
##  Min.   :-1.6491  
##  1st Qu.:-0.8449  
##  Median :-0.1509  
##  Mean   : 0.0000  
##  3rd Qu.: 0.7701  
##  Max.   : 2.5147

summary(x)

##       age               sex               bmi               map         
##  Min.   :-2.2517   Min.   :-0.9375   Min.   :-1.8958   Min.   :-2.3604  
##  1st Qu.:-0.7833   1st Qu.:-0.9375   1st Qu.:-0.7188   1st Qu.:-0.7698  
##  Median : 0.1130   Median :-0.9375   Median :-0.1530   Median :-0.1191  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.7996   3rd Qu.: 1.0643   3rd Qu.: 0.6562   3rd Qu.: 0.7485  
##  Max.   : 2.3253   Max.   : 1.0643   Max.   : 3.5817   Max.   : 2.7729  
##        tc                ldl               hdl               tch          
##  Min.   :-2.66239   Min.   :-2.4279   Min.   :-2.1484   Min.   :-1.60428  
##  1st Qu.:-0.71920   1st Qu.:-0.6375   1st Qu.:-0.7375   1st Qu.:-0.82936  
##  Median :-0.09074   Median :-0.0802   Median :-0.1383   Median :-0.05444  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.59552   3rd Qu.: 0.6267   3rd Qu.: 0.6155   3rd Qu.: 0.72049  
##  Max.   : 3.23219   Max.   : 4.1745   Max.   : 3.8048   Max.   : 3.88992  
##       ltg               glu          
##  Min.   :-2.6480   Min.   :-2.89311  
##  1st Qu.:-0.6982   1st Qu.:-0.69676  
##  Median :-0.0409   Median :-0.02263  
##  Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.6811   3rd Qu.: 0.58626  
##  Max.   : 2.8056   Max.   : 2.84785

Multivariate Regression

Wir benutzen ein multiples lineares Regressionsmodell:

\[ \begin{aligned} y_i &\sim {\rm N}(\mathbf{x}_i\boldsymbol\beta,\sigma^2)\\ \sigma^2&\sim IG(a,b)\\ {\boldsymbol\beta} &\sim{\rm N}(\mu_0,\boldsymbol\Lambda^{-1}_0) \end{aligned} \]

Welche Priori-Information haben wir über \(\beta=(\beta_1,\ldots,\beta_p)\)? Bisher hatten wir keine:

\[p(\beta_j) \propto \text{const.} \]

Dies können wir als uneigentliche Normalverteilung (hier konjugierte Verteilung) mit Varianz unendlich interpretieren.

Ridge-Regression

Nun haben wir aber viele Kovariablen (in diesem Beispiel noch moderat, in anderen Beispielen können es Tausende sein).
Einige Kovariablen sind vielleicht selbst korreliert. Sie erklären dann den selben Einfluß. Wir können dann Kovariablen streichen, ohne dass das Modell schlechter wird.
Wir nehmen also a priori an, dass einige Kovariablen keinen Einfluß haben: Die \(\beta\) Parameter sollen gleich oder nahe Null sein.
Also nehmen wir a priori an

\[ \beta_j \sim N(0,\tau^{-1}) \quad \forall j \]

Um so größer \(\tau\), umso stärker “drückt” diese Priori die \(\beta\) Richtung Null (den Priori-Erwartungswert).
Ein großes \(\tau\) führt also dazu, dass auch a posteriori viele \(\beta\) gleich Null sind.

Wir wollen \(\tau\) nicht festlegen, sondern wieder mitschätzen.
Dazu brauchen wir eine Priori, z.B.: \(\tau\sim\)Ga\((a,b)\).
Man nennt das Modell Bayesianische Ridge-Regression.
Die Kovariableneffekte werden regularisiert, also durch die Priori in eine bestimmte Richtung gedrückt.
Die Kovariablen wurden zuvor standardisiert, damit \(\tau\) auf alle gleich wirkt (sonst würden Kovariablen mit größererer Varianz stärker regularisiert).

Penalisierte Likelihood

Ridge-Regression ist auch aus dem penalisierten Likelihood-Ansatz bekannt.
Die penalisierte Log-Likelihood lautet:

\[ l_{pen}(\beta)=l(\beta)-\frac{\lambda}{2} \sum_{j=1}^p \beta_j^2 \]

wobei \(l(\beta)\) die Log-Likelihood ist.

Die penalisierte Log-Likelihood entspricht der Log-Full-Conditional

Posteriori

Die Posteriori lautet dann:

\[ \begin{aligned} p(\beta,\tau,\sigma^2|y) &\propto \sigma^{-n}\exp\left(-\frac{1}{2\sigma^2}\sum_i(y_i-\sum_j\beta_jx_{ij})^2\right)\\ &\cdot \tau^{p/2}\exp\left(-\frac{\tau}{2}\sum_j\beta_j^2\right)\\ & \cdot \tau^{a-1}\exp(-\tau b)\\ & \cdot \sigma^2{-a_0-1}\exp(-b_0/\sigma^2)\\ \end{aligned} \]

Damit gilt:

\(\beta|\tau,\sigma^2 \sim N(\hat{\beta},\Sigma)\) mit \(\hat{\beta}=(X'X+\tau I)^{-1}X'y\) und \(\Sigma=(X'X+\tau I)^{-1})\).
\(\tau|\beta \sim Ga(a+p/2,b+\sum\beta_j^2/2)\)
\(\sigma^2|\beta,y \sim IG(a_0+n/2,b+\sum(\epsilon_i^2))\) mit \(\epsilon_i=y_i-\sum_j\beta_jx_{ij}\)

Wir können also einen Gibbs-Sampler verwenden.

MCMC

n<-dim(x)[1]
p<-dim(x)[2]
beta<-rep(0,p)
XX <- t(x)%*%x
Xy <- t(x)%*%y
tau <- 1
sigma2 <- 1
a0 <- 1
b0 <- 0.001
a <- 1
b <- 0.1

beta.save<-array(NA,c(p,500))
tau.save<-rep(NA,500)
sigma2.save<-rep(NA,500)

for (i in 1:1000)
{
  Sigma <- solve(XX+tau*diag(p))
  mu <- Sigma%*%Xy
  beta <- mnormt::rmnorm(1,mu,Sigma)
  
  tau <- rgamma(1, a+p/2, b+sum(beta^2)/2)

  sigma2 <- 1/rgamma(1, a0+n/2, b0+sum((y-x%*%beta)^2))
  
  if (i>500)
  {
    beta.save[,i-500]=beta
    tau.save[i-500]=tau
    sigma2.save[i-500]=sigma2
  }
}

Plot

Folgender Plot zeigt die Mediane und 95%-Kredibilitätsintervalle für die zehn Kovariablen:

beta.qu<-apply(beta.save,1,quantile,probs=c(.025,.5,.975))
plot(beta.qu[2,],pch=19,ylim=range(beta.qu),ylab="beta")
for (i in 1:10)
lines(rep(i,2),beta.qu[c(1,3),i],lwd=2)

Bei Ridge werden die Parameter Richtung Null gedrückt.
Aber: Die Parameter werden nicht genau gleich Null.

Relevante Kovariablen sind diejenigen, bei denen das Kredibilitätsintervall die Null nicht überdeckt:

relevant <- c(which(beta.qu[1,]>0),which(beta.qu[3,]<0))
colnames(x)[relevant]

## [1] "bmi" "map" "ltg" "sex"

summary(t(beta.save[relevant,]))

##        V1               V2                V3                V4          
##  Min.   :0.1224   Min.   :0.03028   Min.   :0.07425   Min.   :-0.30642  
##  1st Qu.:0.2706   1st Qu.:0.15577   1st Qu.:0.24038   1st Qu.:-0.17546  
##  Median :0.3117   Median :0.19384   Median :0.30944   Median :-0.13823  
##  Mean   :0.3101   Mean   :0.19246   Mean   :0.30639   Mean   :-0.13842  
##  3rd Qu.:0.3503   3rd Qu.:0.22964   3rd Qu.:0.36598   3rd Qu.:-0.10323  
##  Max.   :0.4416   Max.   :0.34001   Max.   :0.60348   Max.   : 0.03176

Lasso

Alternative zu Ridge ist der Lasso.
In Penalisierungs-Ansätzen auch \(L_1\)-Penalisierung genannt:

\[ pen(\beta)=\sum_j|\beta_j| \]

Für den Bayesianischen Lasso benutzt man folgende Priori:

\[ p(\beta_j) \propto \exp\left(-\frac{\lambda}{2} \sum |\beta_j| \right) \]

Das entspricht einer Laplace-Verteilung mit Erwartungswert 0:

Die Laplace-Verteilung ist sehr nahe der Null höher, drückt also noch mehr genau zum Punkt Null.
Andererseits sind die tails der Dichte höher, also ab etwa \(\beta>2\) ist die Dichte höher. Das heißt, dass a priori höhere \(\beta\)-Werte wahrscheinlicher sind als bei der Ridge-Regression.
Insgesamt regularisiert die Laplace-Priori kleinere \(\beta\)-Wert stärker, aber größere \(\beta\)-Werte schwächer.

Bayesianischer Lasso

Aber: (erstmal) ist kein Gibbs-Sampler mehr möglich: Die Normalverteilungspriori im Ridge war semi-konjugiert zur Normalverteilung der \(y\).
Nach Park, Trevor and Casella, George. The Bayesian Lasso. Journal of American Statistical Association. 103(482):681-686. 2008 gilt aber für:

\[ \begin{aligned} \beta_j |\tau_j^2 &\sim N(0,\tau_j^2)\\ \tau_j|\lambda^2 & \sim Exp(\lambda^2/2)\\ \lambda^2 &\sim Ga(a,b) \end{aligned} \]

dass \(p(\beta|\lambda)=\int p(\beta|\tau,\lambda) d\tau_1,\ldots,d\tau_p)\) genau der gewünschten Laplace-Dichte entspricht.

Auxiliary Augmentation

Solche Ansätze werden allgemein als Auxiliary Augmentation (übersetzt etwa: Hilfserweiterungen) bezeichnet.
Dabei werden zusätzliche (auxiliary) Parameter ins Modell aufgenommen, ohne dass sich die marginale Priori (und damit die marginale Posteriori) ändert.
Diese Erweiterungen (augmentation) sind so gewählt, dass ein Gibbs-Sampler benutzt werden kann.

Implementation

gibbsBLasso = function(x, y, max.steps = 100000) {
  n <- nrow(x)
  m <- ncol(x)
  
  XtX <- t(x) %*% x #Time saving
  xy <- t(x) %*% y
  
  r <- 1
  delta <- 0.1 # 1.78
  
  betaSamples <- matrix(0, max.steps, m)
  sigma2Samples <- rep(0, max.steps)
  invTau2Samples <- matrix(0, max.steps, m)
  lambdaSamples <- rep(0, max.steps)
  
  beta <- drop(backsolve(XtX + diag(nrow=m), xy))
  residue <- drop(y - x %*% beta)
  sigma2 <- drop((t(residue) %*% residue) / n)
  invTau2 <- 1 / (beta * beta)
  lambda <- m * sqrt(sigma2) / sum(abs(beta))
  
  k <- 0
  while (k < max.steps) {
    k <- k + 1
    
    if (k %% 1000 == 0) {
      cat('Iteration:', k, "\r")
    }
    
    # sample beta
    invD <- diag(invTau2)
    invA <- solve(XtX + invD)
    mean <- invA %*% xy
    varcov <- sigma2 * invA
    beta <- drop(mnormt::rmnorm(1, mean, varcov))
    betaSamples[k,] <- beta
    
    # sample sigma2
    shape <- (n+m-1)/2
    residue <- drop(y - x %*% beta)
    scale <- (t(residue) %*% residue + t(beta) %*% invD %*% beta)/2
    sigma2 <- 1/rgamma(1, shape, 1/scale)
    sigma2Samples[k] <- sigma2
    
    # sample tau2
    muPrime <- sqrt(lambda^2 * sigma2 / beta^2)
    lambdaPrime <- lambda^2
    invTau2 <- rep(0, m)
    for (i in seq(m)) {
      invTau2[i] <- VGAM::rinv.gaussian(1, muPrime[i], lambdaPrime)
    }
    invTau2Samples[k, ] <- invTau2
    
    # update lambda
    shape = r + m/2
    scale = delta + sum(1/invTau2)/2
    lambda <- rgamma(1, shape, 1/scale)
    # if (k %% 10 == 0) {
    # low <- k - 9
    # high <- k
    # lambda <- sqrt( 2*m / sum(colMeans(invTau2Samples[low:high, ])) )
    # }
    lambdaSamples[k] <- lambda
  }
  
  betaQu<-apply(betaSamples[seq(max.steps/2, max.steps, 5), ],2,quantile,c(.025,.5,.975))
  return(betaQu)
}

Ausführung

beta.L<-gibbsBLasso(x, y, max.steps = 25000)

## Iteration: 1000 
Iteration: 2000 
Iteration: 3000 
Iteration: 4000 
Iteration: 5000 
Iteration: 6000 
Iteration: 7000 
Iteration: 8000 
Iteration: 9000 
Iteration: 10000 
Iteration: 11000 
Iteration: 12000 
Iteration: 13000 
Iteration: 14000 
Iteration: 15000 
Iteration: 16000 
Iteration: 17000 
Iteration: 18000 
Iteration: 19000 
Iteration: 20000 
Iteration: 21000 
Iteration: 22000 
Iteration: 23000 
Iteration: 24000 
Iteration: 25000

plot(beta.L[2,],pch=19,ylim=range(beta.L),ylab="beta")
for (i in 1:p)
lines(rep(i,2),beta.L[c(1,3),i],lwd=2)

Die Kredibilitätsintervalle der \(\beta\) sind hier sehr viel kleiner
Nur die Kovariablen 1 und 7 stellen sich hier als irrelevant heraus.
Auch hier wird kein \(\beta\) genau gleich Null.

Indikatorvariablen

Problem bei den bisherigen Ansätzen war, dass \(\beta_j\) ein stetiger Parameter ist. Daher gilt automatisch \(p(\beta|\cdot)=0\)!
Idee: Wir müssen mit einer gewissen Priori-Wahrscheinlichkeit \(\beta=0\) zulassen

Indikatorvariablen

Priori-Ansatz:

\[ \beta_j = I_j\tilde{\beta}_j \]

wobei \(I_j\) eine 0/1-Indikatorvariable ist.

Ist \(I_j=0\), wird \(\beta_j\) auf \(0\) gesetzt.
Ist \(I_j=1\), ist \(\beta_j=\tilde{\beta}_j\) mit Priori-Verteilung von \(\tilde{\beta_j}\), z.B. Normalverteilung mit großer Varianz.

Spike and Slab

Computational einfacher ist folgende Priori-Idee:

Ist \(I_j=1\), ist \(\beta_j\) normalverteilt mit großer Varianz (Slab).
Ist \(I_j=0\), hat \(\beta_j\) normalverteilt mit Erwartungswert Null und sehr kleiner Varianz (und damit praktisch gleich Null) (Spike).

Man nennt dies Spike and Slab-Priori.

Spike and Slab-Priori-Modell:

\[ \begin{aligned} \beta|\gamma,\tau^2 &\sim N(0,\tau^2\gamma)\\ \gamma|w & \sim wI_1(\gamma) + (1-w)I_{\nu_0}(\gamma)\\ \tau^2&\sim IG(a_\tau, b_\tau)\\ w&\sim Beta(a_w, b_w) \end{aligned} \]

WGRR

Ishwaran und Rao (2014) zeigen, dass die Spike and Slab prior ein Spezialfall der gewichteten generalisierten Ridge-Regression sind (weighted generalized Ridge regression, WGRR)
Beim WGRR können die \(\beta\) Parameter auf Null gesetzt werden, wir erhalten \(p(\beta=0|y)\)
Wegen des Bayesian Model Averaging gehen aber z.B. bei \(E(y^*|y)\) weiterhin viele/alle Regressionsparameter ein

Implementiert ist der entsprechende MCMC-Algorithmus im spikeSlabGAM-Paket:

library("spikeSlabGAM")
sas<-spikeAndSlab(y,cbind(rep(1,442),x),family="gaussian")

## 
## Model has  11  coefficients in  11 model terms.
## Blockwise sampling: alpha: 1  block(s); xi: 1 block(s).

## Using 2 parallel processes.
## Use 'options(mc.cores = <YourNumberHere>)' to override next time.

## 
## starting chain(s):
## bbb0----------------------------100%
## bb....................b..........

beta<-rbind(sas$samples$beta[[1]],sas$samples$beta[[2]],sas$samples$beta[[3]])
beta.sas<-apply(beta,2,quantile,probs=c(.025,.5,.975))
plot(beta.sas[2,2:11],pch=19,ylim=range(beta.sas),ylab="beta")
for (i in 2:11)
lines(rep(i-1,2),beta.sas[c(1,3),i],lwd=2)

relevant <- c(which(beta.sas[1,]>0),which(beta.sas[3,]<0))
colnames(x)[relevant]

## [1] "map" "tc"  "glu" "bmi"

summary(t(beta[relevant,]))

##        V1                 V2                  V3                 V4            
##  Min.   :-0.15244   Min.   :-0.228669   Min.   :-0.15066   Min.   :-0.1147378  
##  1st Qu.:-0.05234   1st Qu.:-0.001412   1st Qu.:-0.05296   1st Qu.:-0.0418967  
##  Median :-0.01408   Median : 0.002503   Median : 0.03355   Median : 0.0007096  
##  Mean   : 0.04106   Mean   : 0.033415   Mean   : 0.04293   Mean   : 0.0692301  
##  3rd Qu.: 0.12969   3rd Qu.: 0.124985   3rd Qu.: 0.08718   3rd Qu.: 0.1780896  
##  Max.   : 0.33720   Max.   : 0.298764   Max.   : 0.38506   Max.   : 0.3703368

Zusammenfassung

Durch die Priori können wir Modelle regularisieren, hier speziell zur Null drücken. Das erlaubt uns, viele Kovariablen mit ins Modell aufzunehmen, ohne “Overfitting” zu erhalten. Gleichzeitig können wir Variablen ausschliessen.

Die drei gezeigten Ansätze sind unterschiedlich in ihren Annahmen, entsprechend erhalten wir hier auch unterschiedliche Ergebnisse.

Variablenselektion

Multivariate Regression

Ridge-Regression

Penalisierte Likelihood

Posteriori

MCMC

Plot

Lasso

Bayesianischer Lasso

Auxiliary Augmentation

Implementation

Ausführung

Indikatorvariablen

Indikatorvariablen

Spike and Slab

WGRR

Zusammenfassung

Weiter

Bayesianische Variablenselektion