**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 = parallel::detectCores())
library(loo)
library(projpred)
SEED=14124869
```

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

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

- 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
- Draw from posterior distribution using Markov Chain Monte Carlo (MCMC).
- Evaluate how well the model fits the data and possibly revise the model.
- 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).

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 with`glm`

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.

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.

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`

.

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