---
title: "HDS 7. Bootstrap and inference for penalized regression"
author: "Matti Pirinen, University of Helsinki"
date: "10.1.2024"
output:
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
set.seed(15)
```

### From prediction to inference

Previously, we have been using the predictive ability in the test data as the criterion
to choose between regression models.
In practice, we have done this via cross-validation that puts the models
into an empirical test against data 
that were not included in the training data set.
This sounds a very reasonable approach to choose a model when the goal is to predict well
unobserved outcomes.

More broadly, we also intuitively expect that the models that are able to best predict the outcome value 
from the input predictors, are also best in "capturing" the underlying phenomenon/system we want to study.
Therefore, we would also like to extract more information about the phenomenon from those 
best models than what the point estimates of the coefficients alone give us.
A single set of point estimates does not, for example, answer questions such as

 * What is our level of confidence that this particular coefficient $\beta_j$ is 
 non-zero (probability of being important), or that its value is within a particular interval $[a,b]$?
 * Are there dependencies between multiple predictors so that
 they seem to be all needed together for the model to work well? 
 And what are such dependencies in numerical terms?
 * Could the predictor $j$ in the chosen model be equally well replaced by another predictor $k$,
 or some combination of several other predictors?
 
Answers to such questions require that we quantify our uncertainty about the
chosen model and about the values of the coefficients within the chosen model. 
In other words, we want to 
do statistical inference in addition to doing predictions.
While there are standard approaches to inference in any one (unpenalized) regression model, 
such as using the variance-covariance matrix of the ML-estimator to derive confidence intervals/regions,
no equally general approach is (yet) available for the penalized models.
Here we will approach the problem from two directions:

 1. Empirical quantification of variability of our estimation method by using data resampling also known as the **bootstrap**.
 2. Estimating the posterior distribution for the model and its parameters.


The first one, the bootstrap, is a very generally applicable idea of using the observed
data set to generate additional possible data sets and quantifying how much variation there is
in the parameter estimates of the method when applied across the generated data sets.
The flavor of bootstrapping is very practical and it can be interpreted both as an
approximate method to estimate the sampling distribution of the parameters of interest 
and as a Bayesian method to estimate the posterior distribution under certain naive 
prior distribution about the observed data.

The second approach, the full posterior distribution of all the unknown quantities, is
a theoretically sound way to quantify uncertainty, but it is often technically complicated
and hence available in special cases rather than in general. A current work on probabilistic
programming (e.g. [Stan](https://mc-stan.org/)) is aiming to bridge this gap and make
the posterior inference more easily available for a wider audience.


## Bootstrap

Let's consider the problem of estimating the correlation coefficient $r$ 
between variables $x$ and $y$
using a sample of size $n=100$. 
We simulate an example where true $r=0.8$.

Let's start by visualizing the **sampling distribution** of the estimator `cor()`
when it is applied to ($R=5000$) independent sets of $n=100$ observations from the population.
This distribution shows what kind of estimates we are likely to observe in this setting.


```{r}
# function to generate a pair of correlated variables from lecture HDS6
gen.cor.pair <- function(r, n){ # r is target correlation, n is sample size
 u = rnorm(n, 0, sqrt(abs(r))) # temporary variable that is used for generating x1 and x2
 x1 = scale(u + rnorm(n, 0, sqrt(1 - abs(r))))
 x2 = scale(sign(r)*u + rnorm(n, 0, sqrt(1 - abs(r)))) 
 cbind(x1, x2) # returns matrix
}

n = 100
r.true = 0.8
R = 5000 # replicates of data set
r.est = replicate(R, cor(gen.cor.pair(r = r.true, n = n))[1,2]) # returns cor(x,y) for R indep. data sets
hist(r.est, breaks = 50, main = paste0("Sampling distribution at r = ",r.true), col = "orange", 
     xlab = expression(hat(r)), prob = T, sub = paste0(
       "mean = ",round(mean(r.est),3),"; sd = ",round(sd(r.est),3),"; 95%I = (", 
       round(quantile(r.est,c(0.025)),3),",",
       round(quantile(r.est,c(0.975)),3),")"))
abline(v = r.true, col = "red", lwd = 1.5)
```

We see that with $n=100$ the estimates vary from 0.60 to 0.90, 
with the mean close to the true value (in red)
but the distribution showing some skew to left. 

In real life, we do not have a luxury of 5000
independent data sets from which to compute a mean and an interval, 
but instead we have only one data set of $n$ observations.
Let's generate such a single data set.

```{r}
X = gen.cor.pair(r = r.true, n = n)
r.est = cor(X)[1,2]
r.est
```

This particular data set gives a point estimate of $\widehat{r}=0.764$. 
By observing this value alone
we have no idea how much the estimate could vary, if we happened to had sampled 
another set of similar size from the population.

If we are familiar with inference about correlation coefficient, we know that we could estimate
a 95% CI using [Fisher's transformation](https://en.wikipedia.org/wiki/Fisher_transformation) as
```{r}
tanh(atanh(r.est) + c(-1, 1) * 1.96 * 1/sqrt(n-3)) #9 5%CI for cor from Fisher's transf. 
```

However, Fisher's transformation is a specific trick applicable only to
the correlation coefficient estimation and next
we will use the bootstrap as a more general approach 
to derive a confidence interval.
This happens as follows

1. Resample $B$ data sets of size $n$, each **with replacement** from the original data set $X$. 
These $B$ resampled data sets are called **bootstrap samples**, 
and denoted by $X^*_1,\ldots,X_B^*$ to distinguish them from the original data set $X$.

2. Use the distribution of the statistic of interest over the bootstrap samples as
a (posterior) distribution for the statistic of interest.

Let's first see how this works in our case, 
and then discuss more about the motivation and assumptions
behind the procedure.

```{r, eval = T}
B = 5000 #bootstrap samples
r.boot = replicate(B, cor(X[sample(1:n, n, replace = TRUE),])[1,2])
hist(r.boot, breaks = 50, main = paste0("Bootstrap distribution at r.est = ", round(r.est,3)), 
     col = "limegreen", xlab = expression(hat(r)), prob = TRUE, 
     sub = paste0(
       "mean = ",round(mean(r.boot),3),
       "; sd = ",round(sd(r.boot),3),
       "; 95%I = (", round(quantile(r.boot,c(0.025)),3),",",round(quantile(r.boot,c(0.975)),3),")"))
abline(v = r.true, col = "red", lwd = 1.5)
abline(v = r.est, col = "gold", lwd = 1.5)
legend("topleft", lwd = 1.5, col = c("red","gold"), leg = c("true r","est r"), lty = c(1,1))
```

By interpreting this distribution as a posterior distribution of $r$
we can say that a 95% *credible interval* is (0.672,0.835).
Often exactly the same interval is interpreted as a *confidence interval* so
both Bayesian and frequentist interpretations are used for this bootstrap interval.

**(Philosophical) Discussion:** Even though both Bayesian and frequentist interpretations 
are used for this interval, there seems not to be good arguments why this interval is 
a suitable confidence interval,
other than acknowledging that it is actually a credible interval that has a good frequentist calibration.
For example, a nice treatment of the bootstrap by [Hesterberg (2014)](https://arxiv.org/abs/1411.5279) p.54 says that "*The bootstrap percentile interval has no particular derivation - it just works. This is uncomfortable for a mathematically-trained statistician, and unsatisfying for a mathematical statistics course.*"
This "problem" may exist with the frequentist interpretation, but does not exist with the 
Bayesian interpretation, because there this percentile interval of the bootstrap sample 
is the most obvious estimate for the credible interval.

We see that in our example case this interval agrees well with the Fisher's transformation 
above (0.667,0.835),
but let's also do some empirical testing to see how this interval is 
calibrated over possible other data sets that we could have
used for making the bootstrap sample.
Let's make 2000 replicates of the whole bootstrap procedure, 
each replicate coming from its own data set,
and compute the proportion of the replicates
where the 95% bootstrap interval contains the true value of $r$.

```{r, eval = T}
B = 1000
R = 2000
ci.boot = matrix(NA, ncol = 2, nrow = R)
for(ii in 1:R){
  X = gen.cor.pair(r = r.true, n = n) # generate one pair of predictors
  r.boot = replicate(B, cor(X[sample(1:n, n, replace = T),])[1,2])
  ci.boot[ii,] = quantile(r.boot,c(0.025,0.975))
}
# Count in which proportion of data sets the bootstrap 95%CI contained the true r 
mean(ci.boot[,1] <= r.true & r.true <= ci.boot[,2])
```

We have pretty well-calibrated 95% intervals (in frequentist sense)
from the bootstrap samples, meaning that ~95% of the intervals indeed contain the true
parameter. (To be precise, here the intervals are a bit narrow with coverage of about 93.4%
but it still gives a good idea where 95% of the mass is.) 

*To conclude*: Based on the single data set we estimated $\widehat{r}=0.764$ 
and with an application of the bootstrap we found a 95% credible (or confidence) 
interval of (0.67,0.84).


#### What exactly is the bootstrap used for ?

The term "bootstrap" (likely) originates from the saying 
["to pull oneself up by one's bootstraps"](https://en.wiktionary.org/wiki/pull_oneself_up_by_one%27s_bootstraps), 
which is used when a (nearly) impossible task is needed to be completed. 
In statistics, Bradley Efron introduced the term in his 1979
[paper](https://projecteuclid.org/euclid.aos/1176344552) about the method.
The "impossible-sounding" task here is to estimate the sampling distribution of the estimator
by using only one available estimate. 
The reason that the task is not quite as impossible as it may sound is that
the single available estimate results from a data set of $n$ data points, 
and the variation of the estimator can be approximated by resampling 
additional data sets of the same size as the original data set ($n$), 
using the available data set as an approximation to the (larger) target population.

We use bootstrap to estimate the uncertainty in our statistic of interest. 
Typical cases are
to use SD of the bootstrap sample to estimate SE of the parameter estimate
and to use the $\alpha\cdot 100$% interval of the bootstrap sample to approximate the corresponding
credible/confidence interval for the parameter.
The bootstrap works better, the more symmetric and unbiased the sampling distribution
is, and the bootstrap can fail if the shape of the distribution very strongly depends on the
parameter value, or if the sample size is very small.

Note that we **do not** use bootstrap to try to somehow make our single 
point estimate more accurate.
Such a task indeed remains impossible!
By resampling from the existing data set 
we cannot create more accuracy about the parameter value compared to 
what we already have extracted from the available data set.
The bootstrap simply helps us to approximate how much our point estimate might vary if we had
observed another data set of similar size from the same population (frequentist interpretation),
or to approximate the posterior distribution of the parameter of interest (Bayesian interpretation).

#### Bootstrap resamples **with replacement**

When the original data is $X=(x_1,\ldots,x_n)$, a 
bootstrap sample $X^*=(x_1^*,\ldots,x_n^*)$ is a sample of size $n$ from the original data, 
where each $x_j^*$ has been sampled 
independently from among the $n$ values of $X$, 
with equal probability $\frac{1}{n}$ given for each data point $x_k$ of $X$.
In other words, for any bootstrap sample $X^*$,
$\textrm{Pr}(x_j^* = x_k) = \frac{1}{n}$ for all pairs $j,k\in\{1,\ldots, n\}$.

In particular, it follows that the bootstrap sample can include a single data point $x_k$ many times
and can miss some other value $x_{k'}$ completely. 
Actually, the probability that a bootstrap sample
misses a particular data point $x_k$ is 
$(1-\frac{1}{n})^n \underset{n\to\infty}{\longrightarrow}e^{-1}\approx 0.368.$ 
Thus, in large samples,
only about 63.2% of the data points are included in any one bootstrap sample.

For example, if we have observed data set 
$X=(0.6,1.3,2.0)$ we could make $3^3=27$ bootstrap samples, including, e.g.,
$(0.6,0.6,0.6)$,$(1.3,2.0,1.3)$ and $(2.0, 0.6, 1.3)$.

#### Why does the bootstrap work?

The non-parametric bootstrap, that we have been considering, can be seen as first
building an estimate for the probability distribution $f(x)$ of the data $x$
in the target population, and second
using that estimate of the distribution for inference about 
the statistic of interest, $\theta = \theta(x)$, in the target distribution.

Thus, the bootstrap does inference
based on the assumption that the empirical distribution $\widehat{f}$ of the data points observed in 
the data set $X$ is a good approximation to the complete target distribution of the data points
that could have been sampled / will be sampled in the future.
This empirical distribution simply puts a point mass of $\frac{1}{n}$ on each observed data point
and assumes that no other data points are possible. It is clear that this is a crude
approximation, and depends on the situation and the sample size how good results it gives.

A Bayesian interpretation of bootstrap, first framed by Rubin in [1981](https://projecteuclid.org/euclid.aos/1176345338),
assumes that we are estimating the discrete probability of observing each possible data point.
We put an uninformative prior distribution on the probabilities of all possible 
data points and then we update the probabilities as we observe new data points.
Under this very naive model, the expected value of the 
posterior distribution after observing $n$ data points is
the bootstrap distribution where each observed point has 
a probability of $\frac{1}{n}$, and there is a probability of 0 
for any so far unobserved data point.
It follows that the posterior distribution for any given statistic 
$\theta = \theta(x)$ can then be approximated by the bootstrap method as described above.
A detailed explanation is given by [Rasmus Bååth](http://www.sumsar.net/blog/2015/04/the-non-parametric-bootstrap-as-a-bayesian-model/).


#### Read from ISL "5.2 The Bootstrap"


## Bootstrap with penalized regression

Let's apply the bootstrap to a regression setting.
Let's predict `medv` (median price in a region) from the other variables
available in the `Boston` data set.
We will first standardize all the variables.

The results from LASSO are:
```{r, message = FALSE }
library(MASS) #include Boston data set
library(glmnet)
y = scale(Boston$medv)
X = scale(as.matrix(Boston[,-which(names(Boston) == "medv")]))

# Let's apply the standard LASSO estimation at lambda.1se
cv.lasso = cv.glmnet(x = X, y = y, alpha = 1)
lasso.coef = as.numeric(coef(cv.lasso$glmnet.fit, s = cv.lasso$lambda.1se))
cv.lasso$lambda.1se
data.frame(lasso.coef, row.names = c("intcept", colnames(X)))
```

Here, four coefficients are set to 0. 
How confident we are about these being 0?
Let's see whether the same four
coefficients are 0 also in all/most/some bootstrap samples.
We will simply bootstrap the rows of the `Boston` data set, fit LASSO model
and collect the `lambda` parameter and all the coefficients.
We will do this also for `lambda.min` simultaneously.

```{r, fig.width = 9.5}
n = length(y)
B = 500 # bootstrap samples
boot.coef.1se = matrix(NA, nrow = B, ncol = ncol(X) + 1) # +1 for the intercept
colnames(boot.coef.1se) = c("intcept", colnames(X))
boot.coef.min = boot.coef.1se
boot.lambda = matrix(NA, nrow = B, ncol = 2)
colnames(boot.lambda) = c("1se","min")
for(ii in 1:B){
  boot.ind = sample(1:n, size = n, replace = TRUE) # bootstrap indexes
  XX = X[boot.ind,]
  yy = y[boot.ind]
  boot.cv.lasso = cv.glmnet(x = XX, y = yy, alpha = 1)
  boot.coef.1se[ii,] = as.numeric(coef(boot.cv.lasso$glmnet.fit, s = boot.cv.lasso$lambda.1se))
  boot.coef.min[ii,] = as.numeric(coef(boot.cv.lasso$glmnet.fit, s = boot.cv.lasso$lambda.min))
  boot.lambda[ii,] = c(boot.cv.lasso$lambda.1se, boot.cv.lasso$lambda.min)
}
prob.nonzero.1se = apply(abs(boot.coef.1se) > 1e-6, 2, mean) # proportion where each coef was > 1e-6
prob.nonzero.min = apply(abs(boot.coef.min) > 1e-6, 2, mean) 

barplot(rbind(prob.nonzero.1se[-1],prob.nonzero.min[-1]), beside = TRUE, 
        col = c("forestgreen", "green1"), 
        main = "Bootstrap probability of effect (dark = 1se, light = min)")
```

Looking at the `lambda.1se` values (dark green), 
we see that none of the four predictors are always 0, 
although they are among the ones having the smallest probabilities of being non-zero: 
`age` and `rad` ~10%, `indus` ~20% and `tax` and `zn` ~40%.
On the other hand, `lstat`, `ptratio` and `rm` seem to be always non-zero.

At `lambda.min` we see that only `age` and `indus` are sometimes set to 
near zero (< 1e-6) whereas in most cases all variables have 
coefficients of larger magnitude.

How much the estimated `lambda` has varied across the bootstrapping?
```{r}
summary(boot.lambda)
```

Not very much, and we see that `lambda.min` is often very close to zero meaning
that the cross-validation suggests that there is little need for shrinkage in these data.

We can now also look at the variation in the values of coefficients. 
```{r, fig.width = 9.5}
par(mfrow = c(2, 1))
par(mar = c(1, 3, 1, 1))
boxplot(boot.coef.1se[,-1], col = "forestgreen", main = "", 
        xaxt = "n", ylim = c(-0.6, 0.4), cex = 0.5)
abline(h = 0, lty = 2)
text(1.7, 0.35, "lambda.1se", cex = 1.4, col = "forestgreen")
par(mar = c(2, 3, 0, 1))
boxplot(boot.coef.min[,-1], col = "green1", main = "", ylim = c(-0.6, 0.4), cex = 0.5)
abline( h = 0, lty = 2)
text(1.7, 0.35, "lambda.min", cex = 1.4, col = "green1")
```

In the 1se plot, `lstat` is having the largest magnitude, followed by `rm` and then `ptratio`.
The first two clearly separate from the rest and hence we can conclude that they
will be the ones that explain the most variance (since for the standardized variables the
squared effect size is directly proportional to the variance explained). 

The results are broadly similar in the min plot, but there are some clear differences:

1. At `lambda.min` the shrinkage is very weak and does not set any parameters 
exactly to 0 like in the `lambda.1se` model.
Even `indus` and `age`, whose medians are 0 also with `lambda.min`, do show some variation
around 0 in the `lambda.min` model.

2. Some coefficients, like `zn`, `nox`, `dis`, `rad` and `tax`, show clearly larger magnitude
average effects in `lambda.min` than in `lambda.1se`.

Let's study the second observation in more detail by making 
scatter plots of the estimated coefficients of `zn` vs. `dis`, and `rad` vs `tax`. 

```{r, fig.width = 9}
par(mfrow = c(1, 2))
cols = rep( c("green1", "forestgreen"), each = B) 
x1 = c(boot.coef.min[,"zn"], boot.coef.1se[,"zn"])
x2 = c(boot.coef.min[,"dis"], boot.coef.1se[,"dis"])
plot(x1, x2, xlab = "zn", ylab = "dis", cex = 0.8, pch = 3, col = cols)
legend("topright", leg = c("1se","min"), col = c("forestgreen", "green1"), pch = 3)
x1 = c(boot.coef.min[,"rad"], boot.coef.1se[,"rad"])
x2 = c(boot.coef.min[,"tax"], boot.coef.1se[,"tax"])
plot(x1, x2, xlab = "rad", ylab = "tax", cex = 0.8, pch = 3, col = cols)
```

We see a trend of a negative correlation for both pairs likely meaning that, from the 
point of view of the model fit, 
one can compensate an increase in the magnitude of the first variable by an increase
in the opposite direction of the second variable.
When the model is more penalized (dark green), then both values shrink toward 0,
but when there is less penalization, then the values can grow together in magnitude
but in opposite directions.

Finally, with the bootstrap, 
we could now report a 95% credible interval from the LASSO model. 
For example, for `lstat` these would be 
```{r}
res = rbind(quantile(boot.coef.1se[,"lstat"], c(0.025, 0.975)),
            quantile(boot.coef.min[,"lstat"], c(0.025, 0.975)))
rownames(res) = c("1se", "min")
res
```

These analyses were simple to do although required some computation. 
However, the bootstrap may be only a crude approximation to the actual posterior distribution, 
and next we will look at more exact results for the posterior distribution of a LASSO model.

## Bayesian lasso with `blasso()` 

In lecture notes 6, we established that the LASSO model, that was derived as the optimum
of a penalized likelihood function, is simultaneously a maximum-a-posteriori (MAP)
estimate of a Bayesian model. The central property of this Bayesian model was that
the prior distribution on all coefficients was the Laplace distribution
that had a sharp peak at 0, which then leads to a situation where the optimum value
maximizing the posterior distribution (or, equivalently, minimizing the penalized likelihood)
has typically set several coefficients exactly to 0. 

Let's now explore the posterior distributions of the coefficients under a **[Bayesian LASSO](https://people.eecs.berkeley.edu/~jordan/courses/260-spring09/other-readings/park-casella.pdf)**
model, formulated by Park & Casella (2008).
The difference from the pure penalized likelihood version of LASSO is that
here we introduce explicit prior distributions for all parameters in the model and use
computational techniques called Markov chain Monte Carlo (MCMC) to estimate the posteriors.
Technically, this is more complicated than the bootstrapping above, but luckily 
Robert Gramacy has implemented it in the [`blasso()`](https://www.rdocumentation.org/packages/monomvn/versions/1.9-10/topics/blasso) function 
of `monomvn` package.


The Bayesian model behind 
[blasso](https://www.rdocumentation.org/packages/monomvn/versions/1.9-10/topics/blasso) 
includes also the parameters $\mu$, $\lambda$ and $\sigma^2$ and hence 
they also need explicit prior distributions.
The Bayesian LASSO is defined hierarchically in such a way that 
first $\mu$, $\lambda$ and $\sigma^2$ are given "uninformative"
prior distributions to implement the idea that for these values the information would hopefully
come almost solely from the data, not from the prior.
Conditional on $\lambda$ and $\sigma^2$, we then define the LASSO prior (Laplace distribution) 
for the coefficients. This prior allows a strong shrinkage to the coefficient compared
to what the data alone favors, and hence this part of the model has a very "informative" prior.
To complete the model, we state that the observed data is a result of Gaussian
errors added to the linear predictor, just like in standard linear model.
In formulas, the model specification is
\begin{eqnarray*}
\mu &\sim& 1 \textrm{ (i.e., a flat prior over the real numbers)}\\
\sigma^2 &\sim& 1/\sigma^2 \textrm{ (i.e., a flat prior for $\log(\sigma^2)$)}\\
\lambda^2 &\sim& \textrm{"uninformative" Gamma distribution}\\
(\beta_j \,|\,\sigma,\lambda) &\sim& \textrm{Laplace}(0, \sigma/\lambda), \textrm{ for all } j=1,\ldots,p. \\
(y_i \, | \, \mu, \pmb{\beta}, \sigma^2, \pmb{X}_i) &\sim& \mathcal{N}(\mu + \pmb{X}_i^T \pmb{\beta}, \sigma^2) \textrm{ for all } i=1,\ldots,n. 
\end{eqnarray*}

With these distributions, the Bayes rule tells that the joint posterior distribution
of all unknown parameters is proportional to the product of prior and likelihood:
$$
\textrm{Pr}(\pmb{\beta},\mu,\sigma^2,\lambda^2\,|\,\pmb{y},\pmb{X}) \propto 
\textrm{Pr}(\pmb{y}\,|\,\pmb{\beta},\mu,\sigma^2,\lambda^2,\pmb{X})
\cdot \textrm{Pr}(\pmb{\beta}\,|\,\sigma^2,\lambda^2)\textrm{Pr}(\mu)\textrm{Pr}(\sigma^2)\textrm{Pr}(\lambda^2),
$$
where the first term on the right hand side is the Gaussian likelihood of the linear model,
and the remaining 4 terms are the prior distributions that are defined hierarchically
(first independent flat priors on $\mu$, $\sigma^2$ and $\lambda^2$, and then the Laplace
prior on $\pmb{\beta}$ conditional on the values of $\sigma^2$ and $\lambda^2$).
`blasso()` uses an MCMC algorithm to numerically sample from this posterior distribution.

Let's run `blasso()`. It requires the number of MCMC iterations `T` as input.
More iterations will lead to more accurate results, but will also take more time.
`RJ = FALSE` gives the Bayesian LASSO and `RJ = TRUE` adds a Bayesian variable selection via
Reversible Jump MCMC. 
For now, we use `RJ = FALSE`. 
MCMC produces one set of parameter values per each iteration 
and these values can be interpreted as samples from the posterior distribution, 
except at the first few iterations,
where the algorithm is still dependent on its initial values.
Therefore, some initial iterations are often discarded as a "burn-in" period of the algorithm. 
Here we discard the first 50 iterations as burn-in.
Let's run `blasso` and 
compare the posteriors to the earlier bootstrap results (`lambda.min`) via boxplots.

```{r, message = FALSE, fig.width = 9.5, fig.height = 6, eval = T}
library("monomvn")
mcmc.niter = 5000 # how many iterations to run?
blasso.fit = blasso(X = X, y = y, T = mcmc.niter, RJ = FALSE, case = "default", verb = 0)
burnin = 50 # this many first iterations will be discarded as an initial "burn-in" period
iters = burnin:mcmc.niter # use these iterations in results
par(mfrow = c(2, 1))
par(mar = c(1, 3, 1, 1))
boxplot(blasso.fit$beta[iters,], col = "dodgerblue", cex = 0.5,
        names = colnames(boot.coef.min)[-1], main = "", xaxt = "n")
abline(h = 0, lty = 2)
text(1.7, 0.4, "Bayesian LASSO", cex = 1.4, col = "dodgerblue")
par(mar = c(2, 3, 0, 1))
boxplot(boot.coef.min[,-1], col = "green1", main = "", cex = 0.5)
abline( h = 0, lty = 2)
text(2, 0.4, "Bootstrap (lambda.min)", cex = 1.4, col = "green1")
```

The results quite similar and we can conclude that the bootstrap
produced a good picture of the posterior distribution of the LASSO model.

The Boston data set didn't provide much of a need for stronger shrinkage, so let's
next see how the Bayesian LASSO works in a more sparse test case.

**Example 7.1** 
Let's generate similar data as in HDS6 with only a couple of important variables
and many null variables. To keep plots neat, let's use only $p=30$ variables of which
first 3 have non-zero effects, each explaining about 5% of variance of the outcome,
and let's consider $n=250$ individuals.

```{r, fig.width = 9.5, eval = T}
set.seed(122) 
p = 30
n = 250 
phi = 0.05 # variance explained by x_1, should be 0 < phi < 1. 
b = rep(c(sqrt(phi / (1-phi)), 0), c(3, p-3)) # effects 1,2,3 are non-zero, see Lect. 0.1 for "phi"
X = matrix(rnorm(n*p), nrow = n) # columns 1,2,3 of X will have effects, other 27 cols are noise 
eps = scale(rnorm(n, 0, 1)) # epsilon, error term
y = scale(X%*%b + eps) # makes y have mean = 0, var = 1

mcmc.niter = 1000 # how many iterations to run?
blasso.fit = blasso(X = X, y = y, T = mcmc.niter, RJ = FALSE, case = "default", verb = 0)
burnin = 50 # this many first iterations will be discarded as an initial "burn-in" period
iters = burnin:mcmc.niter # use these iterations in results
par(mar = c(3, 3, 0.1, 0.1))
boxplot(blasso.fit$beta[iters,], col = "dodgerblue", main = "", xaxt = "n", cex = 0.5)
abline(h = 0, lty = 2)
```

We see that the first three predictors are having larger effects than most others,
but the model does not strongly force the others to stay at 0.
This is because the model has a single shrinkage parameter $\lambda$ that needs to balance
the opposing needs between the first 3 variables and the rest.
An extension of `blasso()` includes variable selection (`RJ = TRUE`). 
That extension includes a prior distribution on the number of
non-zero coefficients, that by default is the uniform distribution on $\{1,\ldots,p\}$,
and hence then the model is explicitly allowing some coefficients
be exactly zero while others may have larger values that may or may not need shrinkage
via $\lambda$ parameter. Let's try this out.

```{r, fig.width = 9.5, eval = T}
blasso.rj.fit = blasso(X = X, y = y, T = mcmc.niter, RJ = TRUE, case = "default", verb = 0)
par(mar = c(3, 3, 0.1, 0.1))
boxplot(blasso.rj.fit$beta[iters,], col = "steelblue", main = "", xaxt = "n", cex = 0.5)
abline(h = 0, lty = 2)
```

Let's plot the posterior median of the absolute values of each variable from with and without 
the RJ-variable selection. Let's add the diagonal and a red line to denote the true effect size
in the simulation. Let's use red for the 3 true effects and gray for the rest.

```{r, fig.width = 4, fig.height = 4, eval = T}
x1 = apply(abs(blasso.rj.fit$beta[iters,]), 2, median)
x2 = apply(abs(blasso.fit$beta[iters,]), 2, median)
par(mar = c(4, 4, 1, 0.3))
plot(x1, x2, xlab = "RJ", ylab = "no RJ", col = rep(c("red", "gray"), c(3, p-3)), 
     pch = 19, xlim =c(0, 0.25), ylim = c(0, 0.25))
legend("topleft", col= c("red", "gray"), pch = 19, leg = c("true eff", "zero eff"))
abline(0, 1, lty = 2)
abline(v = b[1], lty = 1, col = "red") #true value in simulation
```

We see that with the variable selection we estimate better the effects of both the true effects
and the absolute values of the zero effects. This is because as the 
variable selection sets some noise variables to 0, we do not need to
make global $\lambda$ that large to get a good model, 
and hence we need not penalize the true effects as much as we do
without the variable selection. 

### Other methods

Here we took a look at two standard and complementary ways to assess uncertainty in
penalized regression models.

More generally, statistical inference for high dimensional models is very active topic of research.
See, for example, work on [Selective inference](https://github.com/selective-inference/R-software).





