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

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

• 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 Einfluss. Wir können dann Kovariablen streichen, ohne dass das Modell schlechter wird.
• Wir nehmen also a priori an, dass einige Kovariablen keinen Einfluss 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\) nahe 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ößerer 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 \(\Sigma=\left(\sigma^{-2}X'X+\tau I\right)^{-1}\) und \(\hat{\beta}=\sigma^{-2}\,\Sigma X'y\).
• \(\tau|\beta \sim Ga(a+p/2,b+\sum\beta_j^2/2)\)
• \(\sigma^2|\beta,y \sim IG(a_0+n/2,b_0+\sum\epsilon_i^2/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/sigma2+tau*diag(p))
  mu <- Sigma%*%Xy/sigma2
  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)/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.1961   Min.   :0.09556   Min.   :0.1301   Min.   :-0.25573  
##  1st Qu.:0.2881   1st Qu.:0.16832   1st Qu.:0.2828   1st Qu.:-0.16853  
##  Median :0.3156   Median :0.19956   Median :0.3337   Median :-0.14191  
##  Mean   :0.3167   Mean   :0.19732   Mean   :0.3306   Mean   :-0.14252  
##  3rd Qu.:0.3440   3rd Qu.:0.22469   3rd Qu.:0.3737   3rd Qu.:-0.11688  
##  Max.   :0.4418   Max.   :0.29628   Max.   :0.5266   Max.   :-0.04452

• 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, gegeben der Fehlervarianz \(\sigma^2\):

\[ p(\beta|\sigma^2) = \prod_{j=1}^p \frac{\lambda}{2\sigma} \exp\left(-\frac{\lambda}{\sigma} |\beta_j| \right) \] * Das entspricht für jedes \(\beta_j\) 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 |\sigma^2,\tau_j^2 &\sim N(0,\sigma^2\tau_j^2)\\ \tau_j^2|\lambda^2 & \sim Exp(\lambda^2/2) \end{aligned} \] dass \(p(\beta|\sigma^2,\lambda)=\int p(\beta|\sigma^2,\tau^2)\,p(\tau^2|\lambda^2)\, d\tau_1^2\cdots d\tau_p^2\) genau der gewünschten Laplace-Dichte entspricht. _j^2|2 & Exp(^2/2)

• Dazu kommen die Prioris für die restlichen Parameter: die uneigentliche Priori \(p(\sigma^2)\propto 1/\sigma^2\) für die Fehlervarianz und \(\lambda^2\sim Ga(r,\delta)\); in der folgenden Implementation ist \(r=1\) und \(\delta=0.004\).

Data Augmentation

• Solche Ansätze werden allgemein als Data Augmentation (übersetzt Hilfsvariablen) bezeichnet.
• Dabei werden zusätzliche (auxiliary) Parameter ins Modell aufgenommen, ohne dass sich die marginale Priori (und damit die marginale Posteriori) ändert.
• Diese Variablen (augmentation) sind so gewählt, dass ein effizienterer 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.004
  
  betaSamples <- matrix(0, max.steps, m)
  sigma2Samples <- rep(0, max.steps)
  invTau2Samples <- matrix(0, max.steps, m)
  lambdaSamples <- rep(0, max.steps)
  
  beta <- drop(solve(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)
    rate <- (t(residue) %*% residue + t(beta) %*% invD %*% beta)/2
    sigma2 <- 1/rgamma(1, shape, rate) 
    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
    rate = delta + sum(1/invTau2)/2
    lambda <- sqrt(rgamma(1, shape, rate))
    # 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)

• Auch hier wird kein \(\beta\) genau gleich Null.

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

beta.x <- beta[,2:11] # ohne Intercept
colnames(beta.x) <- colnames(x)
relevant <- c(which(beta.sas[1,2:11]>0),which(beta.sas[3,2:11]<0))
colnames(x)[relevant]

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

summary(beta.x[,relevant,drop=FALSE])

##       bmi              map                 ltg               sex          
##  Min.   :0.1727   Min.   :-0.006272   Min.   :-0.5065   Min.   :-0.26315  
##  1st Qu.:0.3064   1st Qu.: 0.159412   1st Qu.: 0.2508   1st Qu.:-0.14642  
##  Median :0.3356   Median : 0.189919   Median : 0.3281   Median :-0.12146  
##  Mean   :0.3365   Mean   : 0.189780   Mean   : 0.3375   Mean   :-0.11849  
##  3rd Qu.:0.3668   3rd Qu.: 0.220035   3rd Qu.: 0.4122   3rd Qu.:-0.09337  
##  Max.   :0.4818   Max.   : 0.367110   Max.   : 1.7144   Max.   : 0.01211

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 ausschließen.

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

Weiter

• Weiter zur Modellwahl bei hierarchischen Modellen
• Zurück zur Übersicht