Chapter 10
Bayesian machine learning

Lectures 19–20 · shrinkage priors from ridge to the horseshoe, prediction versus explanation, more features than observations, and Bayesian trees and networks

Machine learning, in its dominant tradition, is about prediction. Choose a flexible function, find the parameter values that fit the training data best, usually by minimising a loss plus a penalty whose strength is chosen by cross-validation, and use the result to predict new cases. This works very well, and for many problems it is the right approach. But it leaves something out. A prediction without a statement of how uncertain it is forces the user to guess which predictions to trust. When a wrong prediction is costly, as in a clinical decision, a loan approval or a fraud alert, that uncertainty is not a luxury; it is part of what has to be delivered.

This chapter is about prediction with the uncertainty put back. The models are familiar machine-learning shapes: regression with many features, ensembles of decision trees, neural networks. Each is fitted the Bayesian way, with a prior, a likelihood and a posterior, so that every prediction comes with a distribution rather than a single number. The central idea is that the classical regularisers of machine learning are priors in disguise, and that once we see them as priors we can design better ones.

• Lecture 19 is about shrinkage priors. We will try to price houses from 39 features using only 100 sales, watch an unregularised model overfit, see that ridge regression and the lasso are posterior modes under normal and Laplace priors, and build the horseshoe prior, which shrinks noise hard and leaves real signals almost alone.
• Lecture 20 is about practice. We will compare four priors on 2,830 houses the models never saw, choose a small set of predictors by projection, move to a breast-cancer data set with more features than patients, and look beyond linear models at Bayesian additive regression trees, Gaussian processes and Bayesian neural networks.

Lecture 19 — Shrinkage priors

10.1 What is this house worth?

Dean De Cock (2011) assembled the records of every residential sale in Ames, Iowa, from 2006 to 2010: 2,930 houses, each with about 80 variables describing the lot, the building, its condition and the sale. He published them as a teaching alternative to the much older Boston housing data, and they have since become a standard benchmark in machine learning. We use the version in the modeldata package.

The question is a prediction question: given what is known about a house when it is put on the market, what will it sell for? To make the problem one where regularisation matters, we pretend that we have only a small sample to learn from: a random 100 sales. With 39 candidate predictors, that is fewer than three observations per coefficient. The other 2,830 sales are locked away as a test set, and we will use them only to judge the predictions. That is a luxury real studies rarely have, and it lets us check honestly what each method buys.

library(tidyverse)
library(brms)
library(posterior)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

We build 39 predictors from the raw file: sizes (living area, basement, first floor, garage, lot), counts (rooms, bedrooms, bathrooms, fireplaces, kitchens), ages (year built, year of the last remodel), condition and heating quality on ordinal scales, location (longitude and latitude), and indicators such as central air, a paved drive, a hillside lot, a new-build sale or an abnormal sale (a foreclosure or short sale). Two variables, the month and year of sale, are included even though we expect them to matter little; a good method should discover that.

tmp <- tempfile()
download.file("https://raw.githubusercontent.com/tidymodels/modeldata/main/data/ames.rda", tmp, quiet = TRUE)
load(tmp)                                                   # creates `ames`
q5 <- c("Poor", "Fair", "Typical", "Good", "Excellent")
houses <- ames |> transmute(
  log_price = log(Sale_Price),
  living = Gr_Liv_Area, lot = log(Lot_Area), frontage = Lot_Frontage,
  built = Year_Built, remod = Year_Remod_Add, veneer = Mas_Vnr_Area,
  bsmt = Total_Bsmt_SF, bsmt_unf = Bsmt_Unf_SF, first_fl = First_Flr_SF,
  bsmt_bath = Bsmt_Full_Bath + 0.5 * Bsmt_Half_Bath, full_bath = Full_Bath, half_bath = Half_Bath,
  bedrooms = Bedroom_AbvGr, kitchens = Kitchen_AbvGr, rooms = TotRms_AbvGrd,
  fireplaces = Fireplaces, garage_cars = Garage_Cars, garage_area = Garage_Area,
  deck = Wood_Deck_SF, porch = Open_Porch_SF, enclosed = Enclosed_Porch, screen = Screen_Porch,
  misc = as.integer(Misc_Val > 0), month = Mo_Sold, year_sold = Year_Sold,
  lon = Longitude, lat = Latitude,
  condition = as.numeric(Overall_Cond),                     # ordered: very poor ... very excellent
  central_air = as.integer(Central_Air == "Y"),
  heating_qc = as.numeric(factor(Heating_QC, q5)),
  paved = as.integer(Paved_Drive == "Paved"),
  one_family = as.integer(Bldg_Type == "OneFam"),
  new_sale = as.integer(Sale_Type == "New"),
  abnormal = as.integer(Sale_Condition == "Abnorml"),
  cul_de_sac = as.integer(Lot_Config == "CulDSac"),
  hillside = as.integer(Land_Contour == "HLS"),
  pconc = as.integer(Foundation == "PConc"),                # poured-concrete foundation
  exter_cond = as.numeric(factor(Exter_Cond, q5)),
  func_typ = as.integer(Functional == "Typ"))
c(houses = nrow(houses), predictors = ncol(houses) - 1)

    houses predictors 
      2930         39 

Now the split. One detail matters more than it looks: we standardise every predictor (and the outcome) using the means and standard deviations of the training houses only. Using the whole data set would let information from the test houses leak into the training step, a small version of the leakage of Chapter 9. Standardising also puts every coefficient on the same scale, "change in standardised log price per standard deviation of the predictor", which is what makes a single shrinkage prior for all of them sensible.

set.seed(2010)
idx <- sample(nrow(houses), 100)                             # the 100 training sales
X <- as.matrix(select(houses, -log_price))
centre <- colMeans(X[idx, ]); spread <- apply(X[idx, ], 2, sd)       # training statistics only
y_mean <- mean(houses$log_price[idx]); y_sd <- sd(houses$log_price[idx])
std <- as_tibble(scale(X, centre, spread)) |> mutate(y = (houses$log_price - y_mean) / y_sd)
train <- std[idx, ]; test <- std[-idx, ]
c(train = nrow(train), test = nrow(test), sd_log_price = round(y_sd, 3))
cr <- cor(X[idx, ]); cr[upper.tri(cr, diag = TRUE)] <- NA
as.data.frame(as.table(cr)) |> drop_na() |> arrange(-abs(Freq)) |> head(5) |>
  mutate(Freq = round(Freq, 2))

       train         test sd_log_price 
     100.000     2830.000        0.375 
         Var1        Var2 Freq
1    first_fl        bsmt 0.89
2 garage_area garage_cars 0.89
3       rooms      living 0.81
4       rooms    bedrooms 0.75
5   full_bath      living 0.63

The training sample has 100 houses and the test set 2,830. The standard deviation of log price among the training houses is 0.375, so predicting every house at the training average would miss by about 0.375 on the log scale, roughly 45% in price; that is the baseline any model must beat. Several predictors are strongly correlated: first-floor area with basement area (0.89, because in a one-storey house they share a footprint), garage area with garage capacity in cars (0.89), and total rooms with living area (0.81) and with bedrooms (0.75). Correlated predictors are where coefficient estimates become unstable, because the data cannot easily tell which of two near-duplicates deserves the credit.

10.2 Prediction is not explanation

Every chapter of this course has started with a causal graph, and this one is no exception, but the graph plays a different role. For a prediction question, we do not need the arrows into the outcome to be causal effects, and we do not need to block back-door paths. A variable that is merely correlated with price, for whatever reason, helps predict it. What the graph still decides is which variables are allowed in: anything measured after the outcome (here, after the sale) is leakage, and anything that would not be available when the prediction is needed is useless. The sale condition and new-build indicators pass this test, because both are known when a house is listed; the final negotiated price of the previous sale of the same house, had it been in the file, would not.

The graph also warns us how not to read the fitted model. In a regression with many predictors, each coefficient is a conditional association: the difference in predicted price between two houses that differ in this predictor and are identical in all 38 others. For bedrooms, holding living area fixed, that is the effect of dividing the same floor space into more, smaller rooms, which buyers may well dislike; it is not the effect of adding a bedroom, which usually adds floor area too. For other predictors the coefficient may be confounded by unmeasured quality $U$. Westreich and Greenland (2013) called the habit of interpreting every coefficient in a multivariable model as a causal effect the Table 2 fallacy, after the table in which epidemiology papers traditionally print them. A good prediction model can be a poor guide to what would happen if we intervened, and that is not a flaw in the model; it is a different question, and Chapter 13 is about how to answer it.

10.3 Too many knobs

Start with the model that has no regularisation at all: a linear regression on all 39 standardised predictors with flat priors on the coefficients. Its posterior mean is essentially the least-squares estimate. We fit it, and then score it twice: on the 100 houses it was fitted to, and on the 2,830 it has never seen. The scores are the root mean squared error of the predicted log price (an error of 0.1 is roughly a 10% error in price) and the average log score, the log of the posterior predictive density at each test house's actual price, which rewards honest uncertainty as well as accurate centres (higher is better).

f <- bf(y ~ .)
base <- c(prior(normal(0, 1), class = Intercept), prior(exponential(1), class = sigma))
m_flat <- brm(f, data = train, prior = base, seed = 10, refresh = 0, file = "fits/ch10_flat")   # flat priors on all 39 slopes
score <- function(m, newdata = test) {          # accuracy on the test houses, in log-price units
  mu <- posterior_epred(m, newdata = newdata); ll <- log_lik(m, newdata = newdata)
  c(rmse = sqrt(mean((newdata$y - colMeans(mu))^2)) * y_sd,
    log_score = mean(matrixStats::colLogSumExps(ll) - log(nrow(ll))) - log(y_sd))
}
rbind(training = score(m_flat, train), test = score(m_flat, test)) |> round(3)
fixef(m_flat)[-1, ] |> as.data.frame() |> arrange(-abs(Estimate)) |> head(8) |> round(2)   # the largest eight

          rmse log_score
training 0.084     0.908
test     0.190     0.360
            Estimate Est.Error  Q2.5 Q97.5
living          0.34      0.09  0.15  0.53
bsmt            0.25      0.09  0.08  0.43
garage_cars     0.15      0.09 -0.01  0.32
built           0.15      0.09 -0.04  0.33
hillside        0.14      0.04  0.07  0.21
condition       0.13      0.06  0.02  0.24
heating_qc      0.12      0.05  0.02  0.23
full_bath       0.10      0.06 -0.01  0.21

On the training houses the flat-prior model looks excellent: a root mean squared error of 0.084, about 8% in price, and a log score of 0.91. On the test houses its error is 0.190, more than twice as large, and its log score falls to 0.36. The largest coefficients are sensible (living area, basement area, garage capacity, year built, a hillside lot, condition, heating quality, full bathrooms), so the damage is not a few absurd estimates. It is spread across all 39: every coefficient is nudged by the noise in these particular 100 sales, and the nudges add up to an error in every prediction.

This is overfitting: with 39 knobs and 100 houses, the model can tune itself to the noise in these particular sales. No coefficient is under any pressure to stay small, and each small error in a coefficient becomes an error in every prediction. The widget lets you watch overfitting happen as predictors are added one at a time.

The training $R^2$ only ever goes up: every new predictor, even pure noise, lets least squares fit the sample a little better. The test $R^2$ rises while real signals are being added and then falls, slowly at first and then steeply as the number of predictors approaches the number of observations, where least squares can fit the training data almost perfectly and predicts new data terribly. Ridge regression degrades far more gently. Increase $n$ and the damage recedes; that is why shrinkage matters most when predictors are many relative to observations, and why our 100-house sample is a good place to see it.

10.4 Ridge and lasso are priors

The classical cures for overfitting add a penalty to the least-squares criterion. Ridge regression (Hoerl and Kennard, 1970) minimises

and the lasso (Tibshirani, 1996) replaces ${\beta }_j^2$ by $|{\beta }_j|$. Both pull coefficients toward zero; the strength $\lambda$ is usually chosen by cross-validation. Now take the Bayesian model $y_i\sim \text{Normal}(\alpha +x_i^{\top }\beta ,\sigma )$ with independent priors ${\beta }_j\sim \text{Normal}(0,\tau )$. Its log posterior is, up to a constant,

Multiply by $-2{\sigma }^2$ and this is exactly the ridge criterion with $\lambda ={\sigma }^2/{\tau }^2$. So the ridge estimate is the posterior mode under a normal prior, and a large penalty is a narrow prior. In the same way, the lasso is the posterior mode under a Laplace (double-exponential) prior, $p({\beta }_j)\propto \exp (-|{\beta }_j|/b)$. Regularisation has always been Bayesian inference with the posterior summarised by its mode. Doing the inference fully adds two things. The first is a whole posterior distribution for every coefficient and every prediction, not just a point. The second is that the penalty strength need not be chosen by cross-validation: give the prior scale $\tau$ its own prior and the data estimate it, exactly as the multilevel models of Chapter 6 estimated the spread of group effects. Ridge with a learned $\tau$ is a multilevel model, in which the "groups" are the coefficients.

The two priors shrink differently, and the clearest way to see it is with one coefficient. Suppose the data alone give an estimate $\widehat{\beta }$ with standard error 1, and plot the posterior mean against $\widehat{\beta }$: the shrinkage profile. A flat prior gives the 45° line (no shrinkage).

The normal prior shrinks every estimate by the same proportion: its profile is a straight line through the origin with slope ${\tau }^2/({\tau }^2+1)$. To tame small noisy estimates it must use a small $\tau$, and then it also crushes large, real effects. The Laplace prior's profile bends: it shrinks small estimates harder than the normal does and large estimates less, and for large estimates it runs roughly parallel to the diagonal. But the gap never closes, because its tails are only exponential: large effects are always pulled in by about the same fixed amount. (The posterior mode under a Laplace prior is exactly zero for small estimates, which is how the lasso selects variables; the posterior mean never is.) The horseshoe, the subject of the next section, does what we actually want: small estimates are shrunk almost to zero, and large estimates are left almost untouched.

10.5 The horseshoe

What prior would express the belief "most of these coefficients are essentially zero, but a few may be large"? It needs a tall spike at zero, to shrink noise hard, and very heavy tails, so that a genuinely large effect is not pulled in. Carvalho, Polson and Scott (2010) built such a prior as a scale mixture of normals:

Each coefficient has its own local scale ${\lambda }_j$, and all share one global scale $\tau$. The global scale pulls everything toward zero; the local scales let individual coefficients escape. The half-Cauchy has so much probability near zero and so much in its far tail that most ${\lambda }_j$ are small (the coefficient is shrunk) while a few can be enormous (the coefficient is left alone). The data decide which coefficient is which.

The name comes from a picture. With standardised predictors, $n$ observations and residual standard deviation $\sigma$, the posterior mean of ${\beta }_j$ given its scale is approximately $(1-{\kappa }_j){\widehat{\beta }}_j$, where

is the shrinkage factor: ${\kappa }_j=0$ means no shrinkage, ${\kappa }_j=1$ means the coefficient is set to zero. When $n{\tau }^2/{\sigma }^2=1$ and ${\lambda }_j$ is half-Cauchy, ${\kappa }_j$ has a $\text{Beta}(\frac{1}{2},\frac{1}{2})$ prior: a U-shaped density piled up at both 0 and 1, shaped like a horseshoe. The prior says that each coefficient is very probably either kept or killed, and only rarely half-shrunk. A normal prior, by contrast, gives every coefficient the same $\kappa$.

The global scale matters a great deal, and the obvious default, $\tau \sim \text{Half-Cauchy}(0,1)$, is a poor one. Piironen and Vehtari (2017) showed that it places most of its prior mass on models in which nearly all coefficients are unshrunk, the opposite of what we want. Their recommendation is to state a guess $p_0$ of the number of relevant predictors and set

which centres the prior for $m_{\text{eff}}$ near $p_0$ while still letting the data move it a long way. The guess is soft: a $p_0$ that is wrong by a factor of two or three changes the results little (Exercise 3). The widget simulates the prior distribution of $m_{\text{eff}}$.

With the default half-Cauchy(0, 1) and our Ames setting, the prior puts most of its mass near $m_{\text{eff}}=p$: it expects almost every coefficient to be large, and the horseshoe barely shrinks. With ${\tau }_0$ from the guess $p_0=5$, the prior for $m_{\text{eff}}$ is centred near five, with a long right tail that lets the data ask for more. Notice how the default's behaviour depends on $n$: the more data, the larger $n{\tau }^2$, and the less any coefficient is shrunk. Making the global scale depend on $n$ and $\sigma$, as ${\tau }_0$ does, is what keeps the meaning of the prior stable.

One refinement completes the picture. The pure horseshoe lets a large coefficient be completely unregularised, which causes trouble when a predictor separates the outcome perfectly (logistic regression) or when the likelihood says little about a coefficient: the posterior then has a very long tail that the sampler struggles with. The regularised horseshoe of Piironen and Vehtari (2017) replaces ${\lambda }_j$ by ${\tilde{\lambda }}_j$ with ${\tilde{\lambda }}_j^2=c^2{\lambda }_j^2/(c^2+{\tau }^2{\lambda }_j^2)$, so that even unshrunk coefficients get a weak $\text{Normal}(0,c)$ "slab" prior. In brms, horseshoe() is the regularised version, with the slab scale $c$ given its own prior, and the argument par_ratio $=p_0/(p-p_0)$ sets ${\tau }_0$ by the rule above.

Lecture 20 — Bayesian machine learning in practice

10.6 Four priors, one hundred houses

We now fit three regularised versions of the 39-predictor model and compare them with the flat-prior fit. The Bayesian ridge and Bayesian lasso use a normal and a Laplace prior with a common scale $s_b$ that is learned from the data, with a half-normal(0, 1) prior; brms has no built-in option for a learned scale on regression coefficients, so we declare the parameter ourselves with stanvar(). The horseshoe uses the guess $p_0=5$, so par_ratio $=5/34$. Horseshoe posteriors have narrow funnels near zero (Chapter 6), so we raise adapt_delta to make the sampler take smaller steps.

sb <- stanvar(scode = "real<lower=0> s_b;", block = "parameters") +          # a learned prior scale
      stanvar(scode = "target += normal_lpdf(s_b | 0, 1);", block = "model")  # half-normal(0, 1)
m_ridge <- brm(f, data = train, prior = c(base, prior(normal(0, s_b), class = b)),
               stanvars = sb, seed = 10, refresh = 0, file = "fits/ch10_ridge")
m_lasso <- brm(f, data = train, prior = c(base, prior(double_exponential(0, s_b), class = b)),
               stanvars = sb, seed = 10, refresh = 0, file = "fits/ch10_lasso")
m_hs <- brm(f, data = train, prior = c(base, prior(horseshoe(par_ratio = 5 / 34), class = b)),
            seed = 10, refresh = 0, control = list(adapt_delta = 0.999), file = "fits/ch10_hs")
rbind(ridge = posterior_summary(m_ridge, variable = "s_b"),
      lasso = posterior_summary(m_lasso, variable = "s_b")) |> round(3)
c(horseshoe_divergences = rstan::get_num_divergent(m_hs$fit))

Warning message:
There were 958 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning message:
Examine the pairs() plot to diagnose sampling problems

    Estimate Est.Error  Q2.5 Q97.5
s_b    0.090     0.012 0.069 0.116
s_b    0.069     0.013 0.049 0.097
horseshoe_divergences 
                    0 

The learned prior scales are small: 0.090 for the ridge (95% interval 0.069 to 0.116) and 0.069 for the lasso (0.049 to 0.097). The data chose the strength of the penalty, and chose it to be strong, without any cross-validation. A note on sampling: with adapt_delta = 0.99 the horseshoe produced 35 divergent transitions, so we raised it to 0.999, which removed them. The sampler then reported that 958 of its 4,000 post-warmup iterations reached the maximum tree depth of 10. That warning is about efficiency, not bias: the sampler needed long trajectories to move between the spike and the tails. The diagnostics that matter are fine (all $\widehat{R}\le 1.005$, smallest bulk effective sample size 918). Both warnings are typical of horseshoe posteriors, and both should be reported.

Now the comparison that matters: how well does each model predict the 2,830 houses it never saw? We also compute the leave-one-out estimate of predictive accuracy (PSIS-LOO, introduced in Chapter 1 and explained in Section 9.7) from the training data alone, which is what we would have to rely on without a test set, together with the number of houses whose Pareto-$\widehat{k}$ exceeds 0.7.

models <- list(flat = m_flat, ridge = m_ridge, lasso = m_lasso, horseshoe = m_hs)
loos <- map(models, ~ suppressWarnings(loo(.x)))
tibble(prior = names(models),
       test_rmse = map_dbl(models, ~ score(.x)[["rmse"]]),
       test_log_score = map_dbl(models, ~ score(.x)[["log_score"]]),
       elpd_loo = map_dbl(loos, ~ .x$estimates["elpd_loo", "Estimate"]),
       se = map_dbl(loos, ~ .x$estimates["elpd_loo", "SE"]),
       k_above_0.7 = map_int(loos, ~ sum(.x$diagnostics$pareto_k > 0.7))) |>
  mutate(across(where(is.double), ~ round(.x, 3)))

# A tibble: 4 x 6
  prior     test_rmse test_log_score elpd_loo    se k_above_0.7
  <chr>         <dbl>          <dbl>    <dbl> <dbl>       <int>
1 flat          0.19           0.36     -47.6  8.38          21
2 ridge         0.184          0.317    -37.6  8.55          10
3 lasso         0.181          0.358    -37.6  8.51          11
4 horseshoe     0.176          0.418    -38.8  8.51           9

On the test houses the four priors line up as the theory suggests for prediction error: flat 0.190, ridge 0.184, lasso 0.181 and horseshoe 0.176. The gains are modest, a reduction of about 7% from flat to horseshoe, because with 100 sales and strong signal even least squares is not hopeless. The log score separates them more: the horseshoe is clearly best (0.418 against 0.360 for flat). One result we did not expect is that the ridge has the worst log score (0.317) despite a smaller error than the flat model; its predictive distributions are less well calibrated, which Exercise 7 investigates. The leave-one-out estimates, computed from the training houses alone, rank the flat model last (elpd −47.6 against −37.6 to −38.8, with standard errors near 8.5) but cannot separate the other three, and they slightly prefer ridge and lasso. They are also unreliable: between 9 and 21 of the 100 houses have Pareto-$\widehat{k}$ above 0.7. Without the test set we could have said only that some regularisation helps. And the test set is only one split: Exercise 1 repeats the comparison on ten more random samples, in which the regularised models always beat the flat prior but the ridge and the horseshoe are tied on average, so the horseshoe's margin here is partly the luck of this split.

The coefficients tell the same story from the inside. The figure shows the 90% posterior interval of every coefficient under the flat prior, the learned-scale ridge and the horseshoe, sorted by the size of the horseshoe estimate.

coef_tab <- function(m, prior) as_draws_df(m) |> select(starts_with("b_"), -b_Intercept) |>
  pivot_longer(everything(), names_to = "coef") |> group_by(coef) |>
  summarise(mean = mean(value), lo = quantile(value, .05), hi = quantile(value, .95)) |>
  mutate(coef = str_remove(coef, "^b_"), prior = prior)
ct <- bind_rows(coef_tab(m_flat, "flat"), coef_tab(m_ridge, "ridge"), coef_tab(m_hs, "horseshoe"))
ord <- ct |> filter(prior == "horseshoe") |> arrange(abs(mean)) |> pull(coef)
ct <- ct |> mutate(coef = factor(coef, ord), prior = factor(prior, c("flat", "ridge", "horseshoe")))
ggplot(ct, aes(y = coef, colour = prior)) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, position = position_dodge(width = 0.7)) +
  geom_point(aes(x = mean), position = position_dodge(width = 0.7), size = 1.1) +
  scale_colour_manual(values = c(flat = harbour[["muted"]], ridge = harbour[["sky"]], horseshoe = harbour[["rose"]])) +
  labs(x = "standardised coefficient", y = NULL, colour = NULL)

Under the flat prior every interval is wide (an average width of 0.18) and 12 of the 39 exclude zero. The learned-scale ridge pulls everything in uniformly: the intervals are narrower (0.13 on average) and, because they are narrower around estimates that are still away from zero, 15 now exclude zero. The horseshoe behaves differently: 16 of its 39 posterior means are within 0.02 of zero, and the 8 intervals that exclude zero belong to the predictors a house buyer would name first: living area (0.39), basement area (0.23), year built (0.21), garage capacity (0.15), a hillside lot, heating quality, condition and latitude. Notice what happened to the correlated pairs: the horseshoe gave the credit to living area rather than total rooms, and to garage capacity rather than garage area, and shrank the partner toward zero. Which member of a pair gets the credit is itself uncertain, and another sample might choose the other.

How many coefficients did the horseshoe end up using? The shrinkage factors of Section 10.5 can be computed for every posterior draw: brms stores each coefficient's total prior scale $\tau {\tilde{\lambda }}_j$ as sdb_*, so ${\kappa }_j=1/(1+n{\text{sdb}}_j^2/{\sigma }^2)$ and $m_{\text{eff}}={\sum }_j(1-{\kappa }_j)$.

dr <- as_draws_df(m_hs)
scales <- as.matrix(select(dr, starts_with("sdb_")))
kappa <- 1 / (1 + nrow(train) * scales^2 / dr$sigma^2)
m_eff <- rowSums(1 - kappa)
quantile(m_eff, c(0.05, 0.5, 0.95)) |> round(1)
sort(colMeans(1 - kappa), decreasing = TRUE)[1:8] |> setNames(str_remove(names(sort(colMeans(1 - kappa), decreasing = TRUE)[1:8]), "sdb_")) |> round(2)

  5%  50%  95% 
17.4 22.7 27.8 
     living        bsmt       built garage_cars    hillside  heating_qc 
       0.99        0.94        0.94        0.90        0.89        0.87 
  condition         lat 
       0.81        0.77 

The posterior median of $m_{\text{eff}}$ is 22.7 (90% interval 17.4 to 27.8). The horseshoe treated about 23 of the 39 coefficients as substantially unshrunk, far more than the five our prior guessed. This is a finding about the data, not a failure of the prior: house prices depend on many things a little, not on a handful of things a lot, and the prior's guess was soft enough for the data to overrule it (the widget in Section 10.5 overlays this posterior on the prior when you press the Ames button). The least-shrunk coefficients, with posterior means of $1-{\kappa }_j$, were living area (0.99), basement area and year built (0.94), garage capacity (0.90), a hillside lot (0.89), heating quality (0.87), condition (0.81) and latitude (0.77).

10.7 Which predictors matter?

The horseshoe does not select variables: every coefficient keeps a posterior, however concentrated near zero. Often, though, we want a small model, because measuring 39 things is expensive, or because a model with eight predictors can be explained and audited. A tempting shortcut is to keep the predictors whose 90% intervals exclude zero. It is a poor rule: with correlated predictors (three size measures, two garage measures), the posterior may be sure that one of a group matters while being unsure which, so each interval crosses zero even though dropping the whole group would hurt.

Projection predictive variable selection (Goutis and Robert, 1998; Piironen, Paasiniemi and Vehtari, 2020) handles this properly. Treat the full horseshoe model as the reference model, the best description of the data we have. For each candidate subset of predictors, find the submodel whose predictions are closest to the reference model's (in KL divergence), a projection that is computed draw by draw, so the submodel inherits the reference model's uncertainty. Add predictors one at a time, choosing at each step the one that brings the submodel closest, and stop when the submodel predicts about as well as the reference. The projpred package does this for brms fits.

library(projpred)
vs <- cv_varsel(m_hs, method = "forward", nterms_max = 12, validate_search = FALSE, verbose = FALSE)
size <- suggest_size(vs)
perf <- summary(vs)$perf_sub; ref <- loo(m_hs)$estimates["elpd_loo", "Estimate"]
ggplot(perf, aes(size, elpd)) +
  geom_hline(yintercept = ref, colour = harbour[["muted"]], linetype = 2) +
  geom_pointrange(aes(ymin = elpd - elpd.se, ymax = elpd + elpd.se), colour = harbour[["sky"]], size = 0.3) +
  geom_line(colour = harbour[["sky"]]) +
  scale_x_continuous(breaks = 0:12) + labs(x = "number of predictors", y = "elpd (LOO)")
ranking(vs)$fulldata[1:size]
invisible(capture.output(prj <- project(vs, nterms = size, verbose = FALSE)))
pl <- proj_linpred(prj, newdata = test, integrated = TRUE)
c(size = size, submodel_test_rmse = round(sqrt(mean((test$y - pl$pred[1, ])^2)) * y_sd, 3))

[1] "first_fl"    "built"       "living"      "heating_qc"  "garage_cars"
[6] "hillside"    "bsmt"       
              size submodel_test_rmse 
             7.000              0.199 

The forward search adds first-floor area, year built, living area, heating quality, garage capacity, a hillside lot and basement area, and suggest_size chooses seven predictors. Then comes the honest check: on the 2,830 test houses the projected seven-predictor model has an RMSE of 0.199, worse than the full horseshoe (0.176) and even worse than the flat-prior model (0.190). Two things explain this. First, with validate_search = FALSE the search itself is not cross-validated, so the curve is optimistic about small models; validate_search = TRUE repeats the search inside each fold and is the setting to use in a paper, at much greater cost. Second, with 100 houses the standard errors are wide, and the default rule picks the smallest model whose performance is not clearly worse than the reference, which here is too small: the posterior of $m_{\text{eff}}$ was already saying that about 23 predictors carry signal. Selection always buys simplicity with some accuracy. How much accuracy is something only held-out data can tell you, and it should be reported.

10.8 More features than patients

The Ames problem has 39 predictors and 100 houses. In genomics the ratio is reversed. The Cancer Genome Atlas (TCGA) measured the expression of about 20,500 genes by RNA sequencing in more than a thousand breast tumours (Cancer Genome Atlas Network, 2012), and the UCSC Xena browser distributes the processed matrix, log-transformed and gene by gene, together with the clinical records (Goldman et al., 2020). About one breast cancer in seven is HER2-positive: the gene ERBB2, which encodes the HER2 receptor, is amplified, and the tumour responds to drugs that target the receptor. Clinically, HER2 status is measured by staining and by counting gene copies, not by sequencing. Our question is how well a tumour's expression profile predicts its clinical HER2 status, and which genes carry the signal.

dir.create("data/tcga", showWarnings = FALSE, recursive = TRUE)
xena <- "https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.BRCA.sampleMap%2F"
for (f in c("HiSeqV2.gz", "BRCA_clinicalMatrix"))
  if (!file.exists(file.path("data/tcga", f))) download.file(paste0(xena, f), file.path("data/tcga", f))
expr <- read_tsv("data/tcga/HiSeqV2.gz", show_col_types = FALSE)       # genes x samples, log2(RSEM + 1)
clin <- read_tsv("data/tcga/BRCA_clinicalMatrix", show_col_types = FALSE, guess_max = 5000)
X <- t(as.matrix(expr[, -1])); colnames(X) <- make.names(expr$sample)
her2 <- clin |> filter(substr(sampleID, 14, 15) == "01",               # primary tumours
                       HER2_Final_Status_nature2012 %in% c("Positive", "Negative"))
X <- X[intersect(rownames(X), her2$sampleID), ]
y <- as.integer(her2$HER2_Final_Status_nature2012[match(rownames(X), her2$sampleID)] == "Positive")
set.seed(10)
tr <- sample(nrow(X), 250)                                              # 250 training tumours
v <- sort(apply(X[tr, ], 2, var), decreasing = TRUE)                    # label-free filter, training data only
genes <- names(v)[1:1000]
Z <- scale(X[, genes], colMeans(X[tr, genes]), apply(X[tr, genes], 2, sd))
gtrain <- as_tibble(Z[tr, ]) |> mutate(her2 = y[tr])
gtest <- as_tibble(Z[-tr, ]) |> mutate(her2 = y[-tr])
c(genes_measured = ncol(X), tumours = nrow(X), positive = sum(y), train = nrow(gtrain),
  train_pos = sum(gtrain$her2), test = nrow(gtest), test_pos = sum(gtest$her2),
  ERBB2_variance_rank = match("ERBB2", names(v)))
gcor <- abs(cor(Z[tr, ], y[tr]))[, 1]
writeLines(paste0("var C10T=", jsonlite::toJSON(list(n = 250, p = 1000, gcor = round(unname(gcor), 3)),
                                                auto_unbox = TRUE), ";"), "data/ch10_tcga.js")
round(head(sort(gcor, decreasing = TRUE), 6), 3)

     genes_measured             tumours            positive               train 
              20530                 766                 114                 250 
          train_pos                test            test_pos ERBB2_variance_rank 
                 34                 516                  80                2717 
    PNMT    FGFR4  C2orf54 B4GALNT2   ABCC12    NTRK3 
   0.521    0.315    0.289    0.281    0.279    0.275 

Of the 766 primary tumours with a recorded HER2 status, 114 (14.9%) are positive. We hold out 516 tumours as a test set and train on 250, of which only 34 are positive. Twenty thousand genes is more than we need for a teaching example, so we keep the 1,000 genes whose expression varies most across the training tumours. This filter never looks at HER2 status, and it uses only the training tumours, so the test set stays untouched by every choice we make.

With 1,000 features and 250 tumours, least squares is not just unstable but undefined: there are infinitely many coefficient vectors that fit the training data perfectly. Even among pure noise features, some will be strongly correlated with the outcome by chance, and the more features there are, the stronger the best chance correlation. The widget makes this concrete.

With 250 cases, a single noise feature has a correlation with the outcome of about $\pm 1/\sqrt{250}\approx 0.06$, but the largest of 1,000 noise correlations is typically about 0.21 (in 90% of simulations between 0.19 and 0.25), and the largest of 20,000 is about 0.26. A naive analysis that ranks genes by correlation and reports the top one will find something "significant" in pure noise. In the real data one gene, PNMT, stands far above anything noise produces (0.52), and 28 of the 1,000 genes exceed the typical noise maximum, so there is genuine signal; but below about 0.3, where all but two of the real genes lie, real and spurious associations are mixed together. Sorting them out is exactly the job of a prior that expects most features to be irrelevant.

We fit two Bayesian logistic regressions. The first gives every coefficient an independent normal(0, 1) prior, a choice that would be reasonable with five predictors. The second uses a regularised horseshoe with a guess of ten relevant genes out of 1,000. On two cores, the first takes about two minutes and the second about four, including compilation.

pr_g <- prior(normal(0, 1.5), class = Intercept)
m_gn <- brm(her2 ~ ., data = gtrain, family = bernoulli(), seed = 10, refresh = 0,
            prior = c(pr_g, prior(normal(0, 1), class = b)), file = "fits/ch10_tcga_n")
m_ghs <- brm(her2 ~ ., data = gtrain, family = bernoulli(), seed = 10, refresh = 0,
             prior = c(pr_g, prior(horseshoe(par_ratio = 10 / 990), class = b)),
             control = list(adapt_delta = 0.99), file = "fits/ch10_tcga_hs")
gscore <- function(p, y = gtest$her2) c(
  accuracy = mean((p > 0.5) == y),
  log_loss = -mean(y * log(p) + (1 - y) * log(1 - p)),
  auc = mean(outer(p[y == 1], p[y == 0], ">") + 0.5 * outer(p[y == 1], p[y == 0], "==")))
p_gn <- colMeans(posterior_epred(m_gn, newdata = gtest))
p_ghs <- colMeans(posterior_epred(m_ghs, newdata = gtest))
rbind(always_no = gscore(rep(mean(gtrain$her2), nrow(gtest))),
      normal_prior = gscore(p_gn), horseshoe = gscore(p_ghs)) |> round(3)

             accuracy log_loss   auc
always_no       0.845    0.433 0.500
normal_prior    0.905    0.342 0.882
horseshoe       0.884    0.293 0.863

diag <- function(m) { d <- summarise_draws(as_draws_df(m, variable = "^b_", regex = TRUE), "rhat", "ess_bulk")
  c(divergent = sum(subset(nuts_params(m), Parameter == "divergent__")$Value),
    max_rhat = round(max(d$rhat), 3), min_ess_bulk = round(min(d$ess_bulk))) }
rbind(normal_prior = diag(m_gn), horseshoe = diag(m_ghs))
calib <- function(p, y = gtest$her2) tibble(p, y) |>
  mutate(bin = cut(p, c(0, 0.1, 0.3, 0.6, 1), include.lowest = TRUE)) |>
  summarise(n = n(), predicted = mean(p), observed = mean(y), .by = bin) |> arrange(bin)
calib(p_ghs) |> mutate(across(where(is.double), ~ round(.x, 3)))

             divergent max_rhat min_ess_bulk
normal_prior         0    1.007         3107
horseshoe            0    1.007          702
# A tibble: 4 x 4
  bin           n predicted observed
  <fct>     <int>     <dbl>    <dbl>
1 [0,0.1]     327     0.044    0.046
2 (0.1,0.3]   123     0.173    0.187
3 (0.3,0.6]    39     0.419    0.487
4 (0.6,1]      27     0.811    0.852

The normal prior classifies 467 of the 516 test tumours correctly (90.5%) and the horseshoe 456 (88.4%), against 84.5% for always answering "negative". On the proper scores the order reverses: the horseshoe's log loss is clearly lower (0.293 against 0.342; 0.433 for the base rate), while the AUCs are close (0.863 against 0.882). With 80 positive test tumours, the AUC difference amounts to about 660 of the 34,880 positive–negative pairs ranked the other way round, which is within noise; the difference in log loss is not, because the normal-prior model is overconfident on the tumours it gets wrong. Both samplers ran cleanly (no divergent transitions, $\widehat{R}\le 1.007$). The horseshoe's probabilities are also well calibrated: in the four bins of the table, the predicted rates (0.04, 0.17, 0.42, 0.81) are close to the observed ones (0.05, 0.19, 0.49, 0.85). Where the two models differ most is in what they say about the genes.

coef_tab <- function(m, label) {
  d <- as_draws_df(m, variable = "^b_", regex = TRUE) |> select(starts_with("b_"), -b_Intercept)
  tibble(coef = sub("^b_", "", names(d)), mean = colMeans(d),
         lo = apply(d, 2, quantile, 0.05), hi = apply(d, 2, quantile, 0.95), prior = label) }
gt <- bind_rows(coef_tab(m_gn, "normal(0, 1)"), coef_tab(m_ghs, "horseshoe"))
top <- gt |> filter(prior == "horseshoe") |> slice_max(abs(mean), n = 15) |> arrange(abs(mean)) |> pull(coef)
gt |> filter(coef %in% top) |> mutate(coef = factor(coef, top)) |>
  ggplot(aes(y = coef, colour = prior)) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, position = position_dodge(width = 0.6)) +
  geom_point(aes(x = mean), position = position_dodge(width = 0.6), size = 1.2) +
  scale_colour_manual(values = c(`normal(0, 1)` = harbour[["sky"]], horseshoe = harbour[["rose"]])) +
  labs(x = "coefficient (log-odds per sd of expression)", y = NULL, colour = NULL)

gt |> summarise(mean_width_90 = mean(hi - lo), largest_abs_mean = max(abs(mean)),
                excluding_zero = sum(lo > 0 | hi < 0), .by = prior) |> mutate(across(where(is.double), ~ round(.x, 2)))
gt |> filter(prior == "horseshoe") |> slice_max(abs(mean), n = 5) |> mutate(across(where(is.double), ~ round(.x, 2)))

# A tibble: 2 x 4
  prior        mean_width_90 largest_abs_mean excluding_zero
  <chr>                <dbl>            <dbl>          <int>
1 normal(0, 1)          3.05             1.81              1
2 horseshoe             0.04             1.65              1
# A tibble: 5 x 5
  coef     mean    lo    hi prior    
  <chr>   <dbl> <dbl> <dbl> <chr>    
1 PNMT     1.65  1.14  2.29 horseshoe
2 SLC26A3 -0.12 -0.94  0.01 horseshoe
3 HMGCLL1 -0.1  -0.81  0.01 horseshoe
4 NTRK3   -0.09 -0.82  0.01 horseshoe
5 CHRNA9   0.05 -0.01  0.47 horseshoe

Under the normal prior, the evidence is spread thinly over a thousand genes: the average 90% interval is 3.05 wide on the log-odds scale, and only one interval excludes zero. The model predicts by adding up many small, uncertain contributions. The horseshoe concentrates the evidence: its average interval is 0.04 wide, and a single gene carries the signal, PNMT, with posterior mean 1.65 (90% interval 1.14 to 2.29); no other coefficient has a posterior mean above 0.12 in absolute value. The result is biologically sensible, and instructive. PNMT is not a cancer gene; it encodes an enzyme of adrenaline synthesis. But it sits on chromosome 17 within a few tens of kilobases of ERBB2, inside the smallest region that is amplified in HER2-positive tumours (Staaf et al., 2010), so when ERBB2 is amplified its neighbour is usually amplified with it; in these data its expression tracks HER2 status more closely than any other gene's. Our filter discarded ERBB2, and the horseshoe found the best surviving proxy for it. A coefficient is a predictive association, and here the association is real but the gene is a passenger: a reminder of Section 10.2, on genomic data.

10.9 Beyond linear models

Shrinkage priors work inside a linear model. Three families of Bayesian models go further, capturing nonlinearities and interactions without our specifying them, and still return a posterior over predictions.

Bayesian additive regression trees (BART; Chipman, George and McCulloch, 2010) write the regression function as a sum of $m$ small trees, $f(x)={\sum }_{k=1}^m{T}_k(x)$, typically with $m=200$. The prior keeps each tree shallow (a split at depth $d$ has probability $0.95(1+d{)}^{-2}$) and each leaf value small, so that every tree is a weak learner contributing a small adjustment. The posterior is explored by MCMC: each iteration updates one tree given all the others, proposing to grow, prune or change a split. In spirit BART is gradient boosting with a prior, and it returns predictive intervals and a count of how often each predictor is used in a split. We fit it to the same 100 houses with the dbarts package and compare its 90% predictive intervals with the horseshoe's.

library(dbarts)
set.seed(10)
bt <- bart(as.matrix(select(train, -y)), train$y, as.matrix(select(test, -y)),
           ndpost = 2000, nskip = 1000, verbose = FALSE)
noise <- matrix(rnorm(length(bt$yhat.test)), nrow(bt$yhat.test)) * bt$sigma
pp <- list(horseshoe = posterior_predict(m_hs, newdata = test), BART = bt$yhat.test + noise)
ll_bart <- dnorm(matrix(test$y, nrow(bt$yhat.test), nrow(test), byrow = TRUE), bt$yhat.test, bt$sigma, log = TRUE)
map_dfr(pp, function(d) {
  lo <- apply(d, 2, quantile, 0.05); hi <- apply(d, 2, quantile, 0.95)
  c(rmse = sqrt(mean((test$y - colMeans(d))^2)) * y_sd,
    coverage_90 = mean(test$y > lo & test$y < hi), width_90 = mean(hi - lo) * y_sd)
}, .id = "model") |>
  mutate(log_score = c(score(m_hs)[["log_score"]], mean(matrixStats::colLogSumExps(ll_bart) - log(nrow(ll_bart))) - log(y_sd))) |>
  mutate(across(where(is.double), ~ round(.x, 3)))
use <- colMeans(bt$varcount / rowSums(bt$varcount))
round(sort(use, decreasing = TRUE)[1:6], 3)

# A tibble: 2 x 5
  model      rmse coverage_90 width_90 log_score
  <chr>     <dbl>       <dbl>    <dbl>     <dbl>
1 horseshoe 0.176       0.888    0.439     0.418
2 BART      0.182       0.843    0.427     0.188
     bsmt    living half_bath  first_fl condition  hillside 
    0.038     0.032     0.032     0.032     0.031     0.030 

BART's test RMSE is 0.182, between those of the ridge (0.184) and the lasso (0.181). Its predictive intervals are too narrow: its 90% intervals contain only 84.3% of the test houses, against 88.8% for the horseshoe, and its log score is far lower (0.188 against 0.418). Its variable-use proportions are nearly flat (basement area 0.038; living area, half bathrooms, first-floor area, condition and a hillside lot 0.030 to 0.032; uniform use would be $1/39\approx 0.026$), so with 100 houses it has not learned which predictors matter. This is a fair result, not a defect of BART: flexible models need data to earn their flexibility, and with 100 sales a linear model with a good prior is hard to beat. Exercise 5 gives BART 1,000 houses.

Gaussian processes, met in Chapter 8, are the other main route to flexible Bayesian regression: a prior directly on the function, through a covariance kernel that says how similar the outcomes of two inputs should be. With an automatic relevance determination kernel (one length scale per input), a GP also learns which inputs matter. Its cost grows with the cube of the number of observations, so approximations are needed beyond a few thousand points; brms offers gp() terms with a low-rank approximation.

Bayesian neural networks put a prior on the weights of a neural network and report the posterior predictive distribution, averaging the predictions of many networks drawn from the posterior: $p(y_{\text{new}}\mid x_{\text{new}},y)=\int p(y_{\text{new}}\mid x_{\text{new}},W)p(W\mid y)dW$. Radford Neal (1996) showed that HMC can sample these posteriors for small networks, and that as a one-hidden-layer network becomes infinitely wide, the prior over functions becomes a Gaussian process. For modern networks with millions of weights, exact sampling is out of reach: the posterior is high-dimensional and multimodal, with exact symmetries (swapping two hidden units leaves the function unchanged). Practice relies on approximations: variational inference on the weights (Chapter 9), Laplace approximations around a trained network, stochastic-gradient MCMC, and deep ensembles, a handful of networks trained from different random starts, which behave like a crude sample from the posterior and are often the best-calibrated option available (Lakshminarayanan, Pritzel and Blundell, 2017). The goal is the one this whole course has pursued: a flexible model that says, in calibrated terms, how much it does not know.

10.10 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Ambroise, C. and McLachlan, G. J. (2002). Selection bias in gene extraction on the basis of microarray gene-expression data. Proceedings of the National Academy of Sciences, 99(10), 6562–6566.
2. Cancer Genome Atlas Network (2012). Comprehensive molecular portraits of human breast tumours. Nature, 490(7418), 61–70.
3. Carvalho, C. M., Polson, N. G. and Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2), 465–480.
4. Chipman, H. A., George, E. I. and McCulloch, R. E. (2010). BART: Bayesian additive regression trees. Annals of Applied Statistics, 4(1), 266–298.
5. Collins, G. S., Reitsma, J. B., Altman, D. G. and Moons, K. G. M. (2015). Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD statement. BMJ, 350, g7594.
6. Collins, G. S., Moons, K. G. M., Dhiman, P. et al. (2024). TRIPOD+AI statement: updated guidance for reporting clinical prediction models that use regression or machine learning methods. BMJ, 385, e078378.
7. De Cock, D. (2011). Ames, Iowa: alternative to the Boston housing data as an end of semester regression project. Journal of Statistics Education, 19(3).
8. Duane, S., Kennedy, A. D., Pendleton, B. J. and Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2), 216–222.
9. Gal, Y. and Ghahramani, Z. (2016). Dropout as a Bayesian approximation: representing model uncertainty in deep learning. Proceedings of the 33rd International Conference on Machine Learning, PMLR 48, 1050–1059.
10. George, E. I. and McCulloch, R. E. (1993). Variable selection via Gibbs sampling. Journal of the American Statistical Association, 88(423), 881–889.
11. Goldman, M. J., Craft, B., Hastie, M., Repečka, K., McDade, F., Kamath, A., Banerjee, A., Luo, Y., Rogers, D., Brooks, A. N., Zhu, J. and Haussler, D. (2020). Visualizing and interpreting cancer genomics data via the Xena platform. Nature Biotechnology, 38(6), 675–678.
12. Goutis, C. and Robert, C. P. (1998). Model choice in generalised linear models: a Bayesian approach via Kullback–Leibler projections. Biometrika, 85(1), 29–37.
13. Hill, J. L. (2011). Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 217–240.
14. Hoerl, A. E. and Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67.
15. Lakshminarayanan, B., Pritzel, A. and Blundell, C. (2017). Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in Neural Information Processing Systems, 30.
16. Lindley, D. V. and Smith, A. F. M. (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society, Series B, 34(1), 1–41.
17. MacKay, D. J. C. (1992). A practical Bayesian framework for backpropagation networks. Neural Computation, 4(3), 448–472.
18. Mitchell, T. J. and Beauchamp, J. J. (1988). Bayesian variable selection in linear regression. Journal of the American Statistical Association, 83(404), 1023–1032.
19. Neal, R. M. (1993). Probabilistic Inference Using Markov Chain Monte Carlo Methods. Technical Report CRG-TR-93-1, Department of Computer Science, University of Toronto.
20. Neal, R. M. (1996). Bayesian Learning for Neural Networks. Lecture Notes in Statistics 118. Springer.
21. Park, T. and Casella, G. (2008). The Bayesian lasso. Journal of the American Statistical Association, 103(482), 681–686.
22. Piironen, J., Paasiniemi, M. and Vehtari, A. (2020). Projective inference in high-dimensional problems: prediction and feature selection. Electronic Journal of Statistics, 14(1), 2155–2197.
23. Piironen, J. and Vehtari, A. (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. Electronic Journal of Statistics, 11(2), 5018–5051.
24. Staaf, J., Jönsson, G., Ringnér, M., Vallon-Christersson, J., Grabau, D., Arason, A., Gunnarsson, H., Agnarsson, B. A., Malmström, P.-O., Johannsson, O. T., Loman, N., Barkardottir, R. B. and Borg, Å. (2010). High-resolution genomic and expression analyses of copy number alterations in HER2-amplified breast cancer. Breast Cancer Research, 12(3), R25.
25. Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58(1), 267–288.
26. Westreich, D. and Greenland, S. (2013). The Table 2 fallacy: presenting and interpreting confounder and modifier coefficients. American Journal of Epidemiology, 177(4), 292–298.
27. Zhang, Y. D., Naughton, B. P., Bondell, H. D. and Reich, B. J. (2022). Bayesian regression using a prior on the model fit: the R2-D2 shrinkage prior. Journal of the American Statistical Association, 117(538), 862–874.