candy.Rmd

---
title: "Bayesian variable selection for candy ranking data"
output:
  html_document: default
  html_notebook: default
  pdf_document: default
---

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

This notebook was inspired by Joshua Loftus' two blog posts "Model selection bias invalidates significance tests" http://joshualoftus.com/post/model-selection-bias-invalidates-significance-tests/ and "A conditional approach to inference after model selection"  http://joshualoftus.com/post/conditional-approach-to-inference-after-model-selection/.

In this notebook we illustrate Bayesian inference for model selection, including PSIS-LOO http://link.springer.com/article/10.1007/s11222-016-9696-4 and projection predictive approach http://link.springer.com/article/10.1007/s11222-016-9649-y which makes decision theoretically justified inference after model selection..

Load libraries.
```{r}
library(rstanarm)
options(mc.cores = parallel::detectCores())
library(loo)
library(tidyverse)
library(bayesplot)
library(projpred)
library(fivethirtyeight)
```

We use candy rankings data from fivethirtyeight package (dataset was originally used in this fivethirtyeight story http://fivethirtyeight.com/features/the-ultimate-halloween-candy-power-ranking/).
```{r}
df <- candy_rankings %>% select(-competitorname)
```

We also create a version where winpercent has been replaced with random draws from a normal distribution. This is our "null" data where covariates do not have any predictive information.
```{r}
dfr <- df %>% select(-winpercent)
n <- nrow(dfr)
ry = rnorm(n)
dfr$ry <- ry
```

### Null data

We start first analysing the "null" data set.

The rstanarm package provides stan_glm which accepts same arguments as glm, but makes full Bayesian inference using Stan (Hamiltonian Monte Carlo No-U-Turn-sampling). Doing variable selection we are anyway assuming that some of the variables are not relevant, and thus it is sensible to use priors which assume some of the covariate effects are close to zero. We use regularized horseshoe prior https://projecteuclid.org/euclid.ejs/1513306866 which has lot of prior mass near 0, but also thick tails allowing relevant effects to not shrunk. 
```{r}
fitrhs <- stan_glm(ry ~ ., data = dfr, prior=hs(), seed=150702646, refresh=0)
```
Let's look at the summary:
```{r}
summary(fitrhs)
```
We didn't get divergences, Rhat's are less than 1.1 and n_eff's are useful (see, e.g., http://mc-stan.org/users/documentation/case-studies/rstan_workflow.html).

```{r}
mcmc_areas(as.matrix(fitrhs),prob_outer = .95)
```

All 95% posterior intervals are overlapping 0, regularized horseshoe prior makes the posteriors concentrate near 0, but there is some uncertainty.

We can easily test whether any of the covariates are useful by using cross-validation to compare to a null model,
```{r}
fit0 <- stan_glm(ry ~ 1, data = dfr, seed=150702646, refresh=0)
```

```{r}
(loorhs <- loo(fitrhs))
(loo0 <- loo(fit0))
compare_models(loo0,loorhs)
```

Based on cross-validation covariates together do not contain any useful information, and there is no need to continue with variable selection. This step of checking whether full mode has any predictive power is often ignored especially when non-Bayesian methods are used. If loo (or AIC as Joshua Loftus demonstrated) would be used for stepwise variable selection it is possible that selection process over a large number of models overfits to the data.

To illustrate the robustness of projpred, we make the projective predictive variable selection using the previous model for "null" data. A fast leave-one-out cross-validation approach http://link.springer.com/article/10.1007/s11222-016-9696-4 is used to choose the model size.
```{r, results='hide'}
fitrhs_cv <- cv_varsel(fitrhs, method='forward', cv_method='LOO', n_loo=n)
```

```{r}
fitrhs_cv$varsel$vind
```

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

And we get a loo-cv based recommendation for the model size to choose
```{r}
(nv <- suggest_size(fitrhs_cv,alpha=0.1))
```
We see that projpred agrees that no variables have useful information.

Next we form the projected posterior for the chosen model.
```{r}
projrhs <- project(fitrhs_cv, nv = nv, ns = 4000)
round(colMeans(as.matrix(projrhs)),1)
round(posterior_interval(as.matrix(projrhs)),1)
```
This looks good as the true values for "null" data are intercept=0, sigma=1.

### Original data

Next we repeat the above analysis with original target variable winpercent.

```{r}
fitrhs <- stan_glm(winpercent ~ ., data = df, prior=hs(), seed=150702646, refresh=0)
```
Let's look at the summary:
```{r}
summary(fitrhs)
```

We didn't get divergences, Rhat's are less than 1.1 and n_eff's are useful (see, e.g., http://mc-stan.org/users/documentation/case-studies/rstan_workflow.html).

```{r}
mcmc_areas(as.matrix(fitrhs),prob_outer = .95)
```

95% posterior interval for chocolateTRUE is not overlapping 0, so maybe there is something useful here.

In case of collinear variables it is possible that marginal posteriors overlap 0, but the covariates can still useful for prediction. With many variables it will be difficult to analyse joint posterior to see which variables are jointly relevant. We can easily test whether any of the covariates are useful by using cross-validation to compare to a null model,
```{r}
fit0 <- stan_glm(winpercent ~ 1, data = df, seed=150702646, refresh=0)
```

```{r}
(loorhs <- loo(fitrhs))
(loo0 <- loo(fit0))
compare_models(loo0,loorhs)
```

Based on cross-validation covariates together do contain useful information. If we need just the predictions we can stop here, but if we want to learn more about the relevance of the covariates we can continue with variable selection.

We make the projective predictive variable selection using the previous model for "null" data. A fast leave-one-out cross-validation approach http://link.springer.com/article/10.1007/s11222-016-9696-4 is used to choose the model size.
```{r, results='hide'}
fitrhs_cv <- cv_varsel(fitrhs, method='forward', cv_method='LOO')
```

```{r}
fitrhs_cv$varsel$vind
```

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

Only one variable seems to be needed to get the same performance as the full model.

And we get a loo-cv based recommendation for the model size to choose
```{r}
(nv <- suggest_size(fitrhs_cv,alpha=0.1))
```
projpred recommends to use just one variable.

Next we form the projected posterior for the chosen model.
```{r}
projrhs <- project(fitrhs_cv, nv = nv, ns = 4000)
round(colMeans(as.matrix(projrhs)),1)
round(posterior_interval(as.matrix(projrhs)),1)
```

```{r}
mcmc_areas(as.matrix(projrhs), 
           pars = c('(Intercept)', names(fitrhs_cv$varsel$vind[1])))
```

In our loo and projpred analysis, we find the chocolateTRUE to have predictive information. Other variables may have predictive power, too, but conditionally on chocolateTRUE other variables do not provide additional information.

<br />


### Appendix: Session information

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

<br />


### Appendix: Licenses

* Code &copy; 2017, Aki Vehtari, licensed under BSD-3.
* Text &copy; 2017, Aki Vehtari, licensed under CC-BY-NC 4.0.