Bayesian Logistic Regression with rstanarm

Setup

Load packages

library(tidyverse)
library(caret)
library(GGally)
library(ggplot2)
library(corrplot)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(rstanarm)
options(mc.cores = 1)
library(loo)
library(projpred)
SEED=14124869

1 Introduction

The introduction to Bayesian logistic regression and rstanarm is from a CRAN vignette by Jonah Gabry and Ben Goodrich. CRAN vignette was modified to this notebook by Aki Vehtari. Instead of wells data in CRAN vignette, Pima Indians data is used. The end of this notebook differs significantly from the CRAN vignette. You can read more about how to use rstanarm in several vignettes at CRAN. This vignette includes also demonstrations of calibration plots and projection predictive variable selection.

Acknowledgements: Preprocessing of Pima Indian dataset is from a noteebok by Lao Zhang

---

2 Bayesian analysis

This vignette explains how to estimate generalized linear models (GLMs) for binary (Bernoulli) and Binomial response variables using the stan_glm function in the rstanarm package.

The four steps of a Bayesian analysis are

1. Specify a joint distribution for the outcome(s) and all the unknowns, which typically takes the form of a marginal prior distribution for the unknowns multiplied by a likelihood for the outcome(s) conditional on the unknowns. This joint distribution is proportional to a posterior distribution of the unknowns conditional on the observed data
2. Draw from posterior distribution using Markov Chain Monte Carlo (MCMC).
3. Evaluate how well the model fits the data and possibly revise the model.
4. Draw from the posterior predictive distribution of the outcome(s) given interesting values of the predictors in order to visualize how a manipulation of a predictor affects (a function of) the outcome(s).

Steps 3 and 4 are covered in more depth by the vignette entitled “How to Use the rstanarm Package”. This vignette focuses on Step 1 when the likelihood is the product of conditionally independent binomial distributions (possibly with only one trial per observation).

2.1 Likelihood

For a binomial GLM the likelihood for one observation \(y\) can be written as a conditionally binomial PMF \[\binom{n}{y} \pi^{y} (1 - \pi)^{n - y},\] where \(n\) is the known number of trials, \(\pi = g^{-1}(\eta)\) is the probability of success and \(\eta = \alpha + \mathbf{x}^\top \boldsymbol{\beta}\) is a linear predictor. For a sample of size \(N\), the likelihood of the entire sample is the product of \(N\) individual likelihood contributions.

Because \(\pi\) is a probability, for a binomial model the link function \(g\) maps between the unit interval (the support of \(\pi\)) and the set of all real numbers \(\mathbb{R}\). When applied to a linear predictor \(\eta\) with values in \(\mathbb{R}\), the inverse link function \(g^{-1}(\eta)\) therefore returns a valid probability between 0 and 1.

The two most common link functions used for binomial GLMs are the logit and probit functions. With the logit (or log-odds) link function \(g(x) = \ln{\left(\frac{x}{1-x}\right)}\), the likelihood for a single observation becomes

\[\binom{n}{y}\left(\text{logit}^{-1}(\eta)\right)^y \left(1 - \text{logit}^{-1}(\eta)\right)^{n-y} = \binom{n}{y} \left(\frac{e^{\eta}}{1 + e^{\eta}}\right)^{y} \left(\frac{1}{1 + e^{\eta}}\right)^{n - y}\]

and the probit link function \(g(x) = \Phi^{-1}(x)\) yields the likelihood

\[\binom{n}{y} \left(\Phi(\eta)\right)^{y} \left(1 - \Phi(\eta)\right)^{n - y},\]

where \(\Phi\) is the CDF of the standard normal distribution. The differences between the logit and probit functions are minor and – if, as rstanarm does by default, the probit is scaled so its slope at the origin matches the logit’s – the two link functions should yield similar results. With stan_glm, binomial models with a logit link function can typically be fit slightly faster than the identical model with a probit link because of how the two models are implemented in Stan. Unless the user has a specific reason to prefer the probit link, we recommend the logit simply because it will be slightly faster and more numerically stable.

In theory, there are infinitely many possible link functions, although in practice only a few are typically used. Other common choices are the cauchit and cloglog functions, which can also be used with stan_glm (every link function compatible withglm will work with stan_glm).

A full Bayesian analysis requires specifying prior distributions \(f(\alpha)\) and \(f(\boldsymbol{\beta})\) for the intercept and vector of regression coefficients. When using stan_glm, these distributions can be set using the prior_intercept and prior arguments. The stan_glm function supports a variety of prior distributions, which are explained in the rstanarm documentation (help(priors, package = 'rstanarm')).

As an example, suppose we have \(K\) predictors and believe — prior to seeing the data — that \(\alpha, \beta_1, \dots, \beta_K\) are as likely to be positive as they are to be negative, but are highly unlikely to be far from zero. These beliefs can be represented by normal distributions with mean zero and a small scale (standard deviation). To give \(\alpha\) and each of the \(\beta\)s this prior (with a scale of 1, say), in the call to stan_glm we would include the arguments prior_intercept = normal(0,1) and prior = normal(0,1).

If, on the other hand, we have less a priori confidence that the parameters will be close to zero then we could use a larger scale for the normal distribution and/or a distribution with heavier tails than the normal like the Student t distribution. Step 1 in the “How to Use the rstanarm Package” vignette discusses one such example.

2.2 Posterior

With independent prior distributions, the joint posterior distribution for \(\alpha\) and \(\boldsymbol{\beta}\) is proportional to the product of the priors and the \(N\) likelihood contributions:

\[f\left(\alpha,\boldsymbol{\beta} | \mathbf{y},\mathbf{X}\right) \propto f\left(\alpha\right) \times \prod_{k=1}^K f\left(\beta_k\right) \times \prod_{i=1}^N { g^{-1}\left(\eta_i\right)^{y_i} \left(1 - g^{-1}\left(\eta_i\right)\right)^{n_i-y_i}}.\]

This is posterior distribution that stan_glm will draw from when using MCMC.

2.3 Logistic Regression

When the logit link function is used the model is often referred to as a logistic regression model (the inverse logit function is the CDF of the standard logistic distribution). As an example, here we will show how to carry out a analysis for Pima Indians data set similar to analysis from Chapter 5.4 of Gelman and Hill (2007) using stan_glm.

3 Diabetes data

# file preview shows a header row
diabetes <- read.csv("diabetes.csv", header = TRUE)

# first look at the data set using summary() and str() to understand what type of data are you working
# with
summary(diabetes)

  Pregnancies        Glucose      BloodPressure    SkinThickness  
 Min.   : 0.000   Min.   :  0.0   Min.   :  0.00   Min.   : 0.00  
 1st Qu.: 1.000   1st Qu.: 99.0   1st Qu.: 62.00   1st Qu.: 0.00  
 Median : 3.000   Median :117.0   Median : 72.00   Median :23.00  
 Mean   : 3.845   Mean   :120.9   Mean   : 69.11   Mean   :20.54  
 3rd Qu.: 6.000   3rd Qu.:140.2   3rd Qu.: 80.00   3rd Qu.:32.00  
 Max.   :17.000   Max.   :199.0   Max.   :122.00   Max.   :99.00  
    Insulin           BMI        DiabetesPedigreeFunction      Age       
 Min.   :  0.0   Min.   : 0.00   Min.   :0.0780           Min.   :21.00  
 1st Qu.:  0.0   1st Qu.:27.30   1st Qu.:0.2437           1st Qu.:24.00  
 Median : 30.5   Median :32.00   Median :0.3725           Median :29.00  
 Mean   : 79.8   Mean   :31.99   Mean   :0.4719           Mean   :33.24  
 3rd Qu.:127.2   3rd Qu.:36.60   3rd Qu.:0.6262           3rd Qu.:41.00  
 Max.   :846.0   Max.   :67.10   Max.   :2.4200           Max.   :81.00  
    Outcome     
 Min.   :0.000  
 1st Qu.:0.000  
 Median :0.000  
 Mean   :0.349  
 3rd Qu.:1.000  
 Max.   :1.000  

str(diabetes)

'data.frame':   768 obs. of  9 variables:
 $ Pregnancies             : int  6 1 8 1 0 5 3 10 2 8 ...
 $ Glucose                 : int  148 85 183 89 137 116 78 115 197 125 ...
 $ BloodPressure           : int  72 66 64 66 40 74 50 0 70 96 ...
 $ SkinThickness           : int  35 29 0 23 35 0 32 0 45 0 ...
 $ Insulin                 : int  0 0 0 94 168 0 88 0 543 0 ...
 $ BMI                     : num  33.6 26.6 23.3 28.1 43.1 25.6 31 35.3 30.5 0 ...
 $ DiabetesPedigreeFunction: num  0.627 0.351 0.672 0.167 2.288 ...
 $ Age                     : int  50 31 32 21 33 30 26 29 53 54 ...
 $ Outcome                 : int  1 0 1 0 1 0 1 0 1 1 ...

Pre-processing

# removing those observation rows with 0 in any of the variables
for (i in 2:6) {
      diabetes <- diabetes[-which(diabetes[, i] == 0), ]
}
# scale the covariates for easier comparison of coefficient posteriors
for (i in 1:8) {
      diabetes[i] <- scale(diabetes[i])
}

# modify the data column names slightly for easier typing
names(diabetes)[7] <- "dpf"
names(diabetes) <- tolower(names(diabetes))

n=dim(diabetes)[1]
p=dim(diabetes)[2]
str(diabetes)

'data.frame':   392 obs. of  9 variables:
 $ pregnancies  : num [1:392, 1] -0.7165 -1.0279 -0.0937 -0.4051 -0.7165 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "Pregnancies"
  ..- attr(*, "scaled:center")= Named num 3.3
  .. ..- attr(*, "names")= chr "Pregnancies"
  ..- attr(*, "scaled:scale")= Named num 3.21
  .. ..- attr(*, "names")= chr "Pregnancies"
 $ glucose      : num [1:392, 1] -1.09 0.466 -1.446 2.41 2.151 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "Glucose"
  ..- attr(*, "scaled:center")= Named num 123
  .. ..- attr(*, "names")= chr "Glucose"
  ..- attr(*, "scaled:scale")= Named num 30.9
  .. ..- attr(*, "names")= chr "Glucose"
 $ bloodpressure: num [1:392, 1] -0.3732 -2.4538 -1.6536 -0.0531 -0.8533 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "BloodPressure"
  ..- attr(*, "scaled:center")= Named num 70.7
  .. ..- attr(*, "names")= chr "BloodPressure"
  ..- attr(*, "scaled:scale")= Named num 12.5
  .. ..- attr(*, "names")= chr "BloodPressure"
 $ skinthickness: num [1:392, 1] -0.584 0.557 0.271 1.508 -0.584 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "SkinThickness"
  ..- attr(*, "scaled:center")= Named num 29.1
  .. ..- attr(*, "names")= chr "SkinThickness"
  ..- attr(*, "scaled:scale")= Named num 10.5
  .. ..- attr(*, "names")= chr "SkinThickness"
 $ insulin      : num [1:392, 1] -0.522 0.101 -0.573 3.256 5.806 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "Insulin"
  ..- attr(*, "scaled:center")= Named num 156
  .. ..- attr(*, "names")= chr "Insulin"
  ..- attr(*, "scaled:scale")= Named num 119
  .. ..- attr(*, "names")= chr "Insulin"
 $ bmi          : num [1:392, 1] -0.71 1.425 -0.297 -0.368 -0.425 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "BMI"
  ..- attr(*, "scaled:center")= Named num 33.1
  .. ..- attr(*, "names")= chr "BMI"
  ..- attr(*, "scaled:scale")= Named num 7.03
  .. ..- attr(*, "names")= chr "BMI"
 $ dpf          : num [1:392, 1] -1.031 5.109 -0.796 -1.057 -0.362 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "DiabetesPedigreeFunction"
  ..- attr(*, "scaled:center")= Named num 0.523
  .. ..- attr(*, "names")= chr "DiabetesPedigreeFunction"
  ..- attr(*, "scaled:scale")= Named num 0.345
  .. ..- attr(*, "names")= chr "DiabetesPedigreeFunction"
 $ age          : num [1:392, 1] -0.967 0.209 -0.477 2.17 2.758 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:392] "4" "5" "7" "9" ...
  .. ..$ : chr "Age"
  ..- attr(*, "scaled:center")= Named num 30.9
  .. ..- attr(*, "names")= chr "Age"
  ..- attr(*, "scaled:scale")= Named num 10.2
  .. ..- attr(*, "names")= chr "Age"
 $ outcome      : int  0 1 1 1 1 1 1 0 1 0 ...

print(paste0("number of observations = ", n))

[1] "number of observations = 392"

print(paste0("number of predictors = ", p))

[1] "number of predictors = 9"

Plot correlation structure

corrplot(cor(diabetes[, c(9,1:8)]))

Make outcome to be factor type and create x and y variables

diabetes$outcome <- factor(diabetes$outcome)
# preparing the inputs
x <- model.matrix(outcome ~ . - 1, data = diabetes)
y <- diabetes$outcome
(model_formula <- formula(paste("outcome ~", paste(names(diabetes)[1:(p-1)], collapse = " + "))))

outcome ~ pregnancies + glucose + bloodpressure + skinthickness + 
    insulin + bmi + dpf + age

4 A Bayesian logistic regression model

A Bayesian logistic regression model can be estimated using the stan_glm function. Here we’ll use a Student t prior with 7 degrees of freedom and a scale of 2.5, which, as discussed above, is a reasonable default prior when coefficients should be close to zero but have some chance of being large.

The formula, data and family arguments to stan_glm are specified in exactly the same way as for glm. We’ve also added the seed (for reproducibility). You can read about other possible arguments in the stan_glm documentation (help(stan_glm, package = 'rstanarm')).

t_prior <- student_t(df = 7, location = 0, scale = 2.5)
post1 <- stan_glm(model_formula, data = diabetes,
                 family = binomial(link = "logit"), 
                 prior = t_prior, prior_intercept = t_prior, QR=TRUE,
                 seed = SEED, refresh=0)

stan_glm returns the posterior distribution for the parameters describing the uncertainty related to unknown parameter values:

pplot<-plot(post1, "areas", prob = 0.95, prob_outer = 1)
pplot+ geom_vline(xintercept = 0)

We can extract corresponding posterior median estimates using ‘coef’ function and to get a sense for the uncertainty in our estimates we can use the posterior_interval function to get Bayesian uncertainty intervals. The uncertainty intervals are computed by finding the relevant quantiles of the draws from the posterior distribution. For example, to compute median and 90% intervals we use:

round(coef(post1), 2)

  (Intercept)   pregnancies       glucose bloodpressure skinthickness 
        -1.01          0.26          1.20         -0.02          0.12 
      insulin           bmi           dpf           age 
        -0.10          0.50          0.40          0.35 

round(posterior_interval(post1, prob = 0.9), 2)

                 5%   95%
(Intercept)   -1.25 -0.78
pregnancies   -0.03  0.55
glucose        0.92  1.51
bloodpressure -0.27  0.24
skinthickness -0.19  0.43
insulin       -0.36  0.16
bmi            0.19  0.82
dpf            0.16  0.64
age            0.04  0.67

For more on posterior_interval and interpreting the parameter estimates from a Bayesian model see Step 2 in the “How to Use the rstanarm Package” vignette.

4.1 Leave-one-out cross-validation

rstanarm supports loo package which implements fast Pareto smoothed leave-one-out cross-validation (PSIS-LOO) (Vehtari, Gelman and Gabry, 2017) to compute expected log predictive density (elpd):

(loo1 <- loo(post1, save_psis = TRUE))

Computed from 4000 by 392 log-likelihood matrix.

         Estimate   SE
elpd_loo   -182.4 11.9
p_loo        11.0  1.4
looic       364.7 23.8
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.9]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

Above we see that PSIS-LOO result is reliable as all Pareto \(k\) estimates are small (k< 0.5) (Vehtari, Gelman and Gabry, 2017, p. @Vehtari+etal:PSIS:2022).

4.2 Comparison to a baseline model

Compute baseline result without covariates.

post0 <- update(post1, formula = outcome ~ 1, QR = FALSE, refresh=0)

Compare to baseline

(loo0 <- loo(post0))

Computed from 4000 by 392 log-likelihood matrix.

         Estimate   SE
elpd_loo   -250.0  6.6
p_loo         1.0  0.0
looic       500.0 13.2
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 0.4]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo_compare(loo0, loo1)

      elpd_diff se_diff
post1   0.0       0.0  
post0 -67.6      11.5  

Covariates contain clearly useful information for predictions.

4.3 Other predictive performance measures

For more easily interpretable predictive performance measures, we next compute posterior predictive probabilities and use them to compute classification error.

# Predicted probabilities
linpred <- posterior_linpred(post1)
preds <- posterior_epred(post1)
pred <- colMeans(preds)
pr <- as.integer(pred >= 0.5)
   
# posterior classification accuracy
round(mean(xor(pr,as.integer(y==0))),2)

[1] 0.78

# posterior balanced classification accuracy
round((mean(xor(pr[y==0]>0.5,as.integer(y[y==0])))+mean(xor(pr[y==1]<0.5,as.integer(y[y==1]))))/2,2)

[1] 0.73

The predictive performance above is overoptimistic. To better estimate the predictive performance for new not yet seen data we next use leave-one-out cross-validation:

# LOO predictive probabilities
ploo=E_loo(preds, loo1$psis_object, type="mean", log_ratios = -log_lik(post1))$value
# LOO classification accuracy
round(mean(xor(ploo>0.5,as.integer(y==0))),2)

[1] 0.78

# LOO balanced classification accuracy
round((mean(xor(ploo[y==0]>0.5,as.integer(y[y==0])))+mean(xor(ploo[y==1]<0.5,as.integer(y[y==1]))))/2,2)

[1] 0.72

In this case it happens that all predicted classes are same as with posterior predictions. We can see the small difference in posterior predictive probabilities and LOO probabilities:

qplot(pred, ploo)

4.4 Calibration of predictions

We can also examine calibration of the posterior and LOO predictive probabilities compared to observed binned event rates. We sort the predicted probabilities in ascending order, bin the corresponding binary observations and plot the mean and uncertainty interval based on binomial distribution.

calPlotData<-calibration(y ~ pred + loopred, 
                         data = data.frame(pred=pred,loopred=ploo,y=y), 
                         cuts=10, class="1")
ggplot(calPlotData, auto.key = list(columns = 2))+
  geom_jitter(data=data.frame(pred=pred,loopred=ploo,y=(as.numeric(y)-1)*100), inherit.aes=FALSE,
              aes(x=loopred*100, y=y), height=2, width=0, alpha=0.3) +
  scale_colour_brewer(palette = "Set1")+
  bayesplot::theme_default(base_family = "sans")

From the above calibration plot we see that the model is otherwise calibrated except that for the observations with highest 10% of the predicted probabilities there are less events than predicted.

The above calibration plot is using independent bins for the observations. The following uses a spline fit for smoothing (with narrower uncertainties).

library(splines)
library(MASS)
ggplot(data = data.frame(pred=pred,loopred=ploo,y=as.numeric(y)-1), aes(x=loopred, y=y)) +
  stat_smooth(method='glm', formula = y ~ ns(x, 5), fullrange=TRUE) +
  geom_abline(linetype = 'dashed') +
  labs(x = "Predicted (LOO)", y = "Observed") +
  geom_jitter(height=0.02, width=0, alpha=0.3) +
  scale_y_continuous(breaks=seq(0,1,by=0.1)) +
  xlim(c(0,1))

Dimitriadis, Gneiting, Jordan (2021) proposed recently a new CORP approach for assessing calibration (or reliability as they call it) that uses nonparametric isotonic regression and the pool-adjacent-violators (PAV) algorithm to estimate conditional event probabilities (CEPs) with automated choice of bins.

library(reliabilitydiag)
rd=reliabilitydiag(EMOS = ploo, y = as.numeric(y)-1)
autoplot(rd)+
  labs(x="Predicted (LOO)",
       y="Conditional event probabilities")+
  bayesplot::theme_default(base_family = "sans", base_size=16)

5 Alternative horseshoe prior on weights.

In this example, with n>>p the difference is small, and thus we don’t expect much difference with a different prior and regularized horseshoe prior (Piironen and Vehtari, 2017) is usually more useful for n<p. 

The global scale parameter for horseshoe prior is chosen as recommended by Piironen and Vehtari (2017). To prepare for projection predictive variable selection we

p0 <- 2 # prior guess for the number of relevant variables
tau0 <- p0/(p-p0) * 1/sqrt(n)
hs_prior <- hs(df=1, global_df=1, global_scale=tau0)
t_prior <- student_t(df = 7, location = 0, scale = 2.5)
post2 <- stan_glm(model_formula, data = diabetes,
                 family = binomial(link = "logit"), 
                 prior = hs_prior, prior_intercept = t_prior,
                 seed = SEED, adapt_delta = 0.999, refresh=0)

We see that the regularized horseshoe prior has shrunk the posterior distribution of irrelevant features closer to zero, without affecting the posterior distribution of the relevant features.

pplot <- plot(post2, "areas", prob = 0.95, prob_outer = 1)
pplot + geom_vline(xintercept = 0)

round(coef(post2), 2)

  (Intercept)   pregnancies       glucose bloodpressure skinthickness 
        -0.97          0.16          1.13          0.00          0.07 
      insulin           bmi           dpf           age 
        -0.01          0.41          0.29          0.33 

round(posterior_interval(post2, prob = 0.9), 2)

                 5%   95%
(Intercept)   -1.21 -0.76
pregnancies   -0.04  0.51
glucose        0.88  1.41
bloodpressure -0.14  0.20
skinthickness -0.07  0.40
insulin       -0.21  0.13
bmi            0.07  0.70
dpf            0.02  0.56
age            0.00  0.64

We compute LOO also for the model with the regularized horseshoe prior. Expected log predictive density is higher, but not significantly. This is not surprising as this is a easy data with n>>p.

(loo2 <- loo(post2))

Computed from 4000 by 392 log-likelihood matrix.

         Estimate   SE
elpd_loo   -181.9 11.1
p_loo         9.6  1.3
looic       363.7 22.3
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.1]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo_compare(loo1, loo2)

      elpd_diff se_diff
post2  0.0       0.0   
post1 -0.5       1.4   

Looking at the pairwise posteriors we can see that, for example, posteriors for age and pregnancies effects are correlating and thus we can’t relay on infering variable relevance by looking at the marginal distributions.

mcmc_pairs(as.matrix(post2), pars = c("pregnancies","age"))

6 Projection predictive variable selection

Next we do variable selection using projection predictive variable selection (Piironen, Paasiniemi and Vehtari, 2020; McLatchie et al., 2023). As the number of observations is large compared to the number of covariates, we estimate the performance using LOO-CV only along the search path (validate_search=FALSE), as we may assume that the overfitting in search is negligible (see more about this in McLatchie et al. (2023)).

varsel2 <- cv_varsel(post2, method='forward', cv_method='loo', validate_search=FALSE)

We can now look at the estimated predictive performance of smaller models compared to the full model.

plot(varsel2, stats = c('elpd', 'pctcorr'), deltas=FALSE, text_angle = 45)

As the estimated predictive performance is not going much above the reference model performance, we know that the use of option validate_search=FALSE was safe (see more in McLatchie et al. (2023)).

We get a LOO based recommendation for the model size and the selected variables.

(nsel<-suggest_size(varsel2))

[1] 3

(vsel<-solution_terms(varsel2)[1:nsel])

[1] "glucose" "age"     "bmi"    

Next we form the projected posterior for the chosen model.

proj2 <- project(varsel2, nv = nsel, ns = 4000)
proj2draws <- as.matrix(proj2)
colnames(proj2draws) <- c("Intercept",vsel)
round(colMeans(proj2draws),1)
round(posterior_interval(proj2draws),1)

mcmc_areas(proj2draws, prob = 0.95, prob_outer = 1,
           pars = c('Intercept', vsel))

The projected posterior can be made predictions in the future (with no need to measure the left out variables).

7 Additive non-linear model

We also test an additive non-linear model using stan_gamm4.

post3 <- stan_gamm4(outcome ~ s(glucose) + s(age) + s(bmi) + s(dpf), data = diabetes,
                 family = binomial(link = "logit"), 
                 seed = SEED, refresh=0)

(loo3 <- loo(post3, save_psis = TRUE))

Computed from 4000 by 392 log-likelihood matrix.

         Estimate   SE
elpd_loo   -175.7 10.8
p_loo         9.1  1.1
looic       351.5 21.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.4]).

All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.

loo_compare(loo1, loo3)

      elpd_diff se_diff
post3  0.0       0.0   
post1 -6.6       3.6   

Based on LOO, non-linear model might be better, but the difference is not big.

Let’s look at the calibration.

gammpreds <- posterior_epred(post3)
gammpred <- colMeans(gammpreds)
loogammpred=E_loo(gammpreds, loo3$psis_object, type="mean", log_ratios = -log_lik(post3))$value
calPlotData<-calibration(y ~ loopred + loogammpred, 
                         data = data.frame(loopred=ploo,loogammpred=loogammpred,y=y),
                         cuts=10, class="1")
ggplot(calPlotData, auto.key = list(columns = 2))

Non-linear model appears to be able to partially model the saturating high probabilities and is thus slighlty better calibrated.

The above calibration plot is using independent bins for the observations. The following uses a spline fit for smoothing (with narrower uncertainties).

ggplot(data = data.frame(loopred=loogammpred,y=as.numeric(y)-1), aes(x=loopred, y=y)) + stat_smooth(method='glm', formula = y ~ ns(x, 5), fullrange=TRUE) + geom_abline(linetype = 'dashed') + ylab(label = "Observed") + xlab(label = "Predicted (LOO)") + geom_jitter(height=0.02, width=0, alpha=0.3) + scale_y_continuous(breaks=seq(0,1,by=0.1)) + xlim(c(0,1))

And the new CORP calibration plot.

rd=reliabilitydiag(EMOS = loogammpred, y = as.numeric(y)-1)
autoplot(rd)+
  labs(x="Predicted (LOO)",
       y="Conditional event probabilities")+
  bayesplot::theme_default(base_family = "sans")

References

McLatchie, Y., Rögnvaldsson, S., Weber, F. and Vehtari, A. (2023) ‘Robust and efficient projection predictive inference’, arXiv preprint arXiv:2306.15581.

Piironen, J., Paasiniemi, M. and Vehtari, A. (2020) ‘Projective inference in high-dimensional problems: Prediction and feature selection’, Electronic Journal of Statistics, 14(1), pp. 2155–2197.

Piironen, J. and Vehtari, A. (2017) ‘Sparsity information and regularization in the horseshoe and other shrinkage priors’, Electronic journal of Statistics, 11(2), pp. 5018–5051. doi: 10.1214/17-EJS1337SI.

Vehtari, A., Gelman, A. and Gabry, J. (2017) ‘Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC’, Statistics and Computing, 27(5), pp. 1413–1432. doi: 10.1007/s11222-016-9696-4.

Vehtari, A., Simpson, D., Gelman, A., Yao, Y. and Gabry, J. (2022) ‘Pareto smoothed importance sampling’, arXiv preprint arXiv:1507.02646. Available at: https://arxiv.org/abs/1507.02646v6.

Licenses

• Code © 2017-2018, Aki Vehtari, licensed under BSD-3.
• Text © 2017-2018, Aki Vehtari, licensed under CC-BY-NC 4.0.
• Introduction Text © 2017, Jonah Gabry and Ben Goodrich source, licensed under GPL 3>

Original Computing Environment

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=en_DK.utf8    
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=fi_FI.utf8    LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Europe/Helsinki
tzcode source: system (glibc)

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] reliabilitydiag_0.2.1 MASS_7.3-60           projpred_2.8.0       
 [4] loo_2.6.0.9000        rstanarm_2.32.1       Rcpp_1.0.12          
 [7] bayesplot_1.11.1      corrplot_0.92         GGally_2.2.1         
[10] caret_6.0-94          lattice_0.22-5        lubridate_1.9.3      
[13] forcats_1.0.0         stringr_1.5.1         dplyr_1.1.4          
[16] purrr_1.0.2           readr_2.1.5           tidyr_1.3.1          
[19] tibble_3.2.1          ggplot2_3.5.0         tidyverse_2.0.0      

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3   tensorA_0.36.2.1     jsonlite_1.8.8      
  [4] magrittr_2.0.3       farver_2.1.1         nloptr_2.0.3        
  [7] rmarkdown_2.25       vctrs_0.6.5          minqa_1.2.6         
 [10] base64enc_0.1-3      htmltools_0.5.7      distributional_0.4.0
 [13] curl_5.2.0           pROC_1.18.5          sass_0.4.8          
 [16] parallelly_1.37.0    StanHeaders_2.32.5   bslib_0.6.1         
 [19] htmlwidgets_1.6.4    plyr_1.8.9           zoo_1.8-12          
 [22] cachem_1.0.8         igraph_2.0.2         mime_0.12           
 [25] lifecycle_1.0.4      iterators_1.0.14     pkgconfig_2.0.3     
 [28] colourpicker_1.3.0   Matrix_1.6-5         R6_2.5.1            
 [31] fastmap_1.1.1        future_1.33.1        shiny_1.8.0         
 [34] digest_0.6.34        colorspace_2.1-0     crosstalk_1.2.1     
 [37] labeling_0.4.3       fansi_1.0.6          timechange_0.3.0    
 [40] mgcv_1.9-1           abind_1.4-5          compiler_4.3.2      
 [43] withr_3.0.0          backports_1.4.1      inline_0.3.19       
 [46] shinystan_2.6.0      ggstats_0.5.1        highr_0.10          
 [49] QuickJSR_1.1.3       pkgbuild_1.4.3       lava_1.7.3          
 [52] gtools_3.9.5         ModelMetrics_1.2.2.2 tools_4.3.2         
 [55] httpuv_1.6.14        future.apply_1.11.1  threejs_0.3.3       
 [58] nnet_7.3-19          glue_1.7.0           nlme_3.1-163        
 [61] promises_1.2.1       grid_4.3.2           checkmate_2.3.1     
 [64] reshape2_1.4.4       generics_0.1.3       recipes_1.0.10      
 [67] gtable_0.3.4         tzdb_0.4.0           class_7.3-22        
 [70] data.table_1.15.0    hms_1.1.3            utf8_1.2.4          
 [73] foreach_1.5.2        pillar_1.9.0         markdown_1.12       
 [76] posterior_1.5.0.9000 later_1.3.2          survival_3.5-8      
 [79] tidyselect_1.2.0     miniUI_0.1.1.1       knitr_1.45          
 [82] gridExtra_2.3        V8_4.4.2             stats4_4.3.2        
 [85] xfun_0.42            hardhat_1.3.1        timeDate_4032.109   
 [88] matrixStats_1.2.0    DT_0.31              rstan_2.32.5        
 [91] stringi_1.8.3        yaml_2.3.8           boot_1.3-28         
 [94] evaluate_0.23        codetools_0.2-19     cli_3.6.2           
 [97] RcppParallel_5.1.7   rpart_4.1.23         shinythemes_1.2.0   
[100] xtable_1.8-4         munsell_0.5.0        jquerylib_0.1.4     
[103] globals_0.16.2       parallel_4.3.2       rstantools_2.4.0    
[106] ellipsis_0.3.2       gower_1.0.1          dygraphs_1.1.1.6    
[109] lme4_1.1-35.1        listenv_0.9.1        ipred_0.9-14        
[112] ggridges_0.5.6       scales_1.3.0         xts_0.13.2          
[115] prodlim_2023.08.28   rlang_1.1.3          shinyjs_2.1.0