# Setup

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 = parallel::detectCores())
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.

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

# 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 "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  "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 "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  "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 "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  "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 "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  "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)]))