mesquite.Rmd

---
title: "Collinear demo with mesquite bushes"
author: "[Aki Vehtari](https://users.aalto.fi/~ave/)"
date: "First version 2018-01-16. 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
link-citations: yes
---

# 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** - Requires projpred 1.0.1+
```{r}
library(rstanarm)
library(arm)
options(mc.cores = parallel::detectCores())
library(loo)
library(ggplot2)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(GGally)
library(projpred)
```

# Introduction

This notebook demonstrates collinearity in multipredictor regression.
Example of predicting the yields of mesquite bushes comes from [Gelman
and Hill (2007)](http://www.stat.columbia.edu/~gelman/arm/).  The
outcome variable is the total weight (in grams) of photosynthetic
material as derived from actual harvesting of the bush. The predictor
variables are:

- diam1: diameter of the canopy (the leafy area of the bush)
in meters, measured along the longer axis of the bush
- diam2: canopy diameter measured along the shorter axis
- canopy height: height of the canopy
- total height: total height of the bush
- density: plant unit density (# of primary stems per plant unit)
- group: group of measurements (0 for the first group, 1 for the second group)

# Data

```{r warning=FALSE,message=FALSE}
dat <- read.table("mesquite.dat", header=TRUE) %>% 
  mutate_if(is.character, as.factor) %>%
  mutate_if(is.factor, as.numeric)
summary(dat)
```

Plot data
```{r}
ggpairs(dat, diag=list(continuous="barDiag"))
```

Additional transformed variables
```{r}
dat$CanVol <- dat$Diam1 * dat$Diam2 * dat$CanHt
dat$CanAre <- dat$Diam1 * dat$Diam2
dat$CanSha <- dat$Diam1 / dat$Diam2
```

It may be reasonable to fit on the logarithmic scale, so that effects are multiplicative rather than additive (we'll return to checking this assumption in another notebook).

# Maxiumum likelihood estimate

We first illustrate the problem with maxiumum likelihood estimate
```{r}
lm1 <- lm(formula = log(LeafWt) ~ log(CanVol) + log(CanAre) + log(CanSha) + log(TotHt) + log(Dens) + Group, data = dat)
display(lm1)
```

GroupMCD seems to be only variable which has coeffiecent far away from zero. Let's try making a model with just the group variable.
```{r}
lm2 <- lm(formula = log(LeafWt) ~ Group, data = dat)
display(lm2)
```
Hmmm.... R-squared dropped a lot, so it seems that other variables are useful even if estimated effects and their standard errors indicate that they are not relevant. There are approach for maximum likelihood estimated models to investigate this, but we'll switch now to Bayesian inference using [`rstanarm`](https://cran.r-project.org/package=rstanarm).

# Bayesian inference

The corresponding `rstanarm` model fit using `stan_glm`
```{r}
fitg <- stan_glm(formula = log(LeafWt) ~ log(CanVol) + log(CanAre) + log(CanSha) + log(TotHt) + log(Dens) + Group, data = dat, refresh=0)
```

Print summary for some diagnostics.
```{r}
summary(fitg)
```

Rhats and n_effs are good (see, e.g., [RStan workflow](http://mc-stan.org/users/documentation/case-studies/rstan_workflow.html)), but QR transformation usually makes sampling work even better (see, [The QR Decomposition For Regression Models](http://mc-stan.org/users/documentation/case-studies/qr_regression.html))

Print summary for some diagnostics.
```{r}
summary(fitg)
```
Use of QR decomposition improved sampling efficiency (actually we get superefficient sampling, ie better than independent sampling) and we continue with this model.

Instead of looking at the tables, it's easier to look at plots
```{r, fig.height=4, fig.width=8}
mcmc_areas(as.matrix(fitg), prob = .5, prob_outer = .95)
```

All 95% posterior intervals except for GroupMCD are overlapping 0 and it seems we have serious collinearity problem.

Looking at the pairwise posteriors we can see high correlations especially between log(CanVol) and log(CanAre).
```{r}
mcmc_pairs(as.matrix(fitg),pars = c("log(CanVol)","log(CanAre)","log(CanSha)","log(TotHt)","log(Dens)"))
```

If look more carefully on of the subplots, we see that although marginal posterior intervals overlap 0, some pairwise joint posteriors are not overlapping 0. Let's look more carefully the joint posterior of log(CanVol) and log(CanAre).
```{r}
mcmc_scatter(as.matrix(fitg), pars = c("log(CanVol)","log(CanAre)")) +
  geom_vline(xintercept=0) +
  geom_hline(yintercept=0)
```

From the joint posterior scatter plot, we can see that 0 is far away fron the typical set.

In case of even more variables with some being relevant and some irrelevant, 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}
fitg0 <- update(fitg, formula = log(LeafWt) ~ 1, QR=FALSE)
```

We compute leave-one-out cross-validation elpd's using PSIS-LOO [@Vehtari+etal:PSIS-LOO:2017]
```{r}
(loog <- loo(fitg))
(loog0 <- loo(fitg0))
```
The model with variables has one bad Pareto $k$ value [@Vehtari+etal:PSIS-LOO:2017]. We can fix that by computing the corresponding leave-one-out-posterior exactly [@Vehtari+etal:PSIS-LOO:2017].
```{r}
(loog <- loo(fitg, k_threshold=0.7))
```

And then we can compare the models.
```{r}
loo_compare(loog0, loog)
```
Based on cross-validation covariates together contain significant information to improve predictions.

We might want to choose some variables 1) because we don't want to observe all the variables in the future (e.g. due to the measurement cost), or 2) we want to most relevant variables which we define here as a minimal set of variables which can provide similar predictions to the full model.

LOO can be used for model selection, but we don't recommend it for variable selection as discussed by @Piironen+Vehtari:2017a. The reason for not using LOO in variable selection is that the selection process uses the data twice, and in case of large number variable combinations the selection process overfits and can produce really bad models. Using the usual posterior inference given the selected variables ignores that the selected variables are conditonal on the selection process and simply setting some variables to 0 ignores the uncertainty related to their relevance.

@Piironen+Vehtari:2017a also show that a projection predictive approach can be used to make a model reduction, that is, choosing a smaller model with some coefficients set to 0. The projection predictive approach solves the problem how to do inference after the selection. The solution is to project the full model posterior to the restricted subspace. See more by @Piironen+etal:projpred:2020.

We make the projective predictive variable selection using the previous full model. A fast leave-one-out cross-validation approach [@Vehtari+etal:PSIS-LOO:2017] is used to choose the model size.
```{r, results='hide'}
fitg_cv <- cv_varsel(fitg, method='forward', cv_method='LOO', nloo = nrow(dat))
```

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

And we get a loo-cv based recommendation for the model size to choose
```{r}
(nsel <- suggest_size(fitg_cv, alpha=0.1))
(vsel <- solution_terms(fitg_cv)[1:nsel])
```
We see that `r nsel` variables is enough to get the same predictive accuracy as with all 4 variables.

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

```{r}
mcmc_areas(projdraws, pars=c("Intercept",vsel), prob_outer=0.99)
```


<br />

# References {.unnumbered}

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

# Licenses {.unnumbered}

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

# Original Computing Environment {.unnumbered}

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

<br />