diabetes.Rmd

---
title: "Bayesian Logistic Regression with rstanarm"
author: "Aki Vehtari, Jonah Gabry, and Ben Goodrich"
date: "First version 2017-07-17. Last modified `r format(Sys.Date())`."
output:
  html_document:
    fig_caption: yes
    toc: TRUE
    toc_depth: 2
    number_sections: TRUE
    toc_float:
      smooth_scroll: FALSE
bibliography: modelsel.bib
csl: harvard-cite-them-right.csl
---

# Setup  {.unnumbered}

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache=FALSE, message=FALSE, error=FALSE, warning=FALSE, comment=NA, out.width='95%')
```

**Load packages**
```{r}
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
```

# Introduction

The introduction to Bayesian logistic regression and **rstanarm** is from a 
[CRAN vignette][1] by Jonah Gabry and Ben Goodrich. [CRAN vignette][1] was modified to this notebook by [Aki Vehtari][2]. *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][3].

Acknowledgements: Preprocessing of Pima Indian dataset is from a [noteebok][4] by [Lao Zhang][5]

  [1]: https://cran.r-project.org/web/packages/rstanarm/vignettes/binomial.html
  [2]: https://users.aalto.fi/~ave/
  [3]: https://cran.r-project.org/web/packages/rstanarm/
  [4]: https://www.kaggle.com/laozhang/d/uciml/pima-indians-diabetes-database/statistical-learning-with-r/run/445129
  [5]: https://www.kaggle.com/laozhang

---

# 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"][2]. This vignette focuses on Step 1 when the likelihood is
the product of conditionally independent binomial distributions (possibly with
only one trial per observation).

## 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](https://en.wikipedia.org/wiki/Logit) and 
[probit](https://en.wikipedia.org/wiki/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.

## 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.
  
## 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)](http://www.stat.columbia.edu/~gelman/arm/) 
using `stan_glm`.

[2]: https://cran.r-project.org/web/packages/rstanarm/vignettes/rstanarm.html


# Diabetes data
```{r}
# 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)
str(diabetes)
```

Pre-processing
```{r}

# 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)
print(paste0("number of observations = ", n))
print(paste0("number of predictors = ", p))
```

Plot correlation structure
```{r}
corrplot(cor(diabetes[, c(9,1:8)]))
```

Make outcome to be factor type and create x and y variables
```{r}
diabetes$outcome <- factor(diabetes$outcome)
# preparing the inputs
x <- model.matrix(outcome ~ . - 1, data = diabetes)
y <- diabetes$outcome
```

# 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')`).

```{r, results='hide'}
t_prior <- student_t(df = 7, location = 0, scale = 2.5)
post1 <- stan_glm(outcome ~ ., data = diabetes,
                 family = binomial(link = "logit"), 
                 prior = t_prior, prior_intercept = t_prior, QR=TRUE,
                 seed = SEED)
```

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

```{r}
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:

```{r}
round(coef(post1), 2)
round(posterior_interval(post1, prob = 0.9), 2)
```

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"][1] vignette.

[1]: https://cran.r-project.org/web/packages/rstanarm/vignettes/rstanarm.html

## Leave-one-out cross-validation

__rstanarm__ supports __loo__ package which implements fast Pareto smoothed leave-one-out cross-validation (PSIS-LOO) [@Vehtari+etal:PSIS-LOO:2017] to compute expected log predictive density (elpd):

```{r}
(loo1 <- loo(post1, save_psis = TRUE))
```

Above we see that PSIS-LOO result is reliable as all Pareto $k$ estimates are small (k< 0.5) [@Vehtari+etal:PSIS-LOO:2017; @Vehtari+etal:PSIS:2019].

## Comparison to a baseline model

Compute baseline result without covariates.
```{r}
post0 <- update(post1, formula = outcome ~ 1, QR = FALSE)
```
Compare to baseline
```{r}
(loo0 <- loo(post0))
```
```{r}
loo_compare(loo0, loo1)
```
Covariates contain clearly useful information for predictions.

## Other predictive performance measures

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

```{r}
# Predicted probabilities
linpred <- posterior_linpred(post1)
preds <- posterior_linpred(post1, transform=TRUE)
pred <- colMeans(preds)
pr <- as.integer(pred >= 0.5)
   
# posterior classification accuracy
round(mean(xor(pr,as.integer(y==0))),2)
# 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)
```

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:

```{r}
# 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)
# 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)
```

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:

```{r}
qplot(pred, ploo)
```

## 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.

```{r}
calPlotData<-calibration(y ~ pred + loopred, 
                         data = data.frame(pred=pred,loopred=ploo,y=as.factor(y)), 
                         cuts=10, class="1")
ggplot(calPlotData, auto.key = list(columns = 2))
```

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).
```{r}
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='gam', formula = y ~ s(x, k=20), method.args = list(family = "binomial")) +
  geom_abline(linetype = 'dashed') +
  labs(x = "Predicted (LOO)", y = "Observed") +
  geom_jitter(aes(x=loopred, y=y*0.96+0.02), height=0.02, width=0, alpha=0.3) +
  scale_y_continuous(breaks=seq(0,1,by=0.1), limits=c(0,1)) +
  xlim(c(0,1))
```

# Alternative horseshoe prior on weights. 

In this example, with n&gt;&gt;p the difference is small, and thus we don't expect much difference with a different prior and regularized horseshoe prior [@Piironen+Vehtari:RHS:2017] is usually more useful for n&lt;p. 

The global scale parameter for horseshoe prior is chosen as recommended by Piironen and Vehtari [-@Piironen+Vehtari:RHS:2017].

```{r, results='hide'}
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(outcome ~ ., data = diabetes,
                 family = binomial(link = "logit"), 
                 prior = hs_prior, prior_intercept = t_prior,
                 seed = SEED, adapt_delta = 0.9)
```

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.

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

```{r}
round(coef(post2), 2)
round(posterior_interval(post2, prob = 0.9), 2)
```

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&gt;&gt;p.

```{r}
(loo2 <- loo(post2))
```

```{r}
loo_compare(loo1, loo2)
```

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.
```{r}
mcmc_pairs(as.matrix(post2), pars = c("pregnancies","age"))
```

# Projection predictive variable selection

Next we do variable selection using projection predictive variable selection [@Piironen+etal:projpred:2018]
```{r, results='hide'}
varsel2 <- cv_varsel(post2, method='forward', nloo = n)
```

```{r}
varsel2$vind
```

We can now look at the estimated predictive performance of smaller models compared to the full model.
```{r}
varsel_plot(varsel2, stats = c('elpd', 'pctcorr'))
```

And we get a LOO based recommendation for the model size to choose
```{r}
(nv <- suggest_size(varsel2, alpha=0.2))
```
Next we form the projected posterior for the chosen model.
```{r}
proj2 <- project(varsel2, nv = nv, ns = 4000)
round(colMeans(as.matrix(proj2)),1)
round(posterior_interval(as.matrix(proj2)),1)
```
```{r}
mcmc_areas(as.matrix(proj2), , prob = 0.95, prob_outer = 1,
           pars = c('(Intercept)', names(varsel2$vind[1:nv])))
```

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

# Additive non-linear model

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

```{r, results='hide'}
post3 <- stan_gamm4(outcome ~ s(glucose) + s(age) + s(bmi), data = diabetes,
                 family = binomial(link = "logit"), 
                 seed = SEED)
```
```{r}
(loo3 <- loo(post3, save_psis = TRUE))
loo_compare(loo1, loo3)
```
Based on LOO, non-linear model might be better, but the difference is not big.

Let's look at the calibration.
```{r, results='hide'}
gammpreds <- posterior_linpred(post3, transform=TRUE)
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).
```{r}
ggplot(data = data.frame(loopred=loogammpred,y=as.numeric(y)-1), aes(x=loopred, y=y)) +
  stat_smooth(method='gam', formula = y ~ s(x, k=20), method.args = list(family = "binomial")) +
  geom_abline(linetype = 'dashed') +
  ylab(label = "Observed") +
  xlab(label = "Predicted (LOO)") +
  geom_jitter(aes(x=loopred, y=y*0.96+0.02), height=0.02, width=0, alpha=0.3) +
  scale_y_continuous(breaks=seq(0,1,by=0.1), limits=c(0,1)) +
  xlim(c(0,1))
```

<br />

# References {.unnumbered}

<div id="refs"></div>

# Licenses {.unnumbered}

* Code &copy; 2017-2018, Aki Vehtari, licensed under BSD-3.
* Text &copy; 2017-2018, Aki Vehtari, licensed under CC-BY-NC 4.0.
* Introduction Text &copy; 2017, Jonah Gabry and Ben Goodrich [source](https://cran.r-project.org/web/packages/rstanarm/vignettes/binomial.html), licensed under GPL 3>

# Original Computing Environment {.unnumbered}

```{r}
sessionInfo()
```

<br />