modelselection.Rmd

---
title: "Model selection"
output:
  pdf_document: default
  html_notebook: default
  html_document:
    df_print: paged
---

In this notebook we're going to analyse different techniques for model selection and afterwards we're going to discuss their shortcomings. 

## Selection criteria

First of all, we're going to look at different criteria to compare models based on their performance and complexity.
```{r}
require(ISLR)
head(Hitters)
summary(Hitters)
```
```{r}
# removing the NA 
dim(Hitters)
Hitters<- na.omit(Hitters)
dim(Hitters)
```
We're going to use cross-validation to compare the results from different selection criteria. 
```{r}
nfolds <- 10 
n <- dim(Hitters)[1]
folds <- cut(1:n, nfolds, labels = F)
# a bit of shuffling
indices <- sample(1:n, size=n, replace=F)
```
```{r}
library(leaps)



get.bss.test.error<- function(train, test, cv.best){
  # estimates the error on the test dataset for the best model 
  # according to each criteria 
  all.best<- regsubsets(x=Salary~.,data=train,nbest=1,
                        nvmax=dim(train)[2]-1, # using all variables
                        method="forward" )
  s <- summary(all.best)
  r2 <- coef(all.best, id=which.max(s$rsq))
  adjr2 <- coef(all.best, id=which.max(s$adjr2))
  cp <- coef(all.best, id=which.min(s$cp))
  bic <- coef(all.best, id=which.min(s$bic))
  cv.coefs <- coef(all.best, id=cv.best)
  # test predictions 
  r2.pred <- model.matrix(Salary~.,test)[,names(r2)]%*%r2
  adjr2.pred <- model.matrix(Salary~.,test)[,names(adjr2)]%*%adjr2
  cp.pred <- model.matrix(Salary~.,test)[,names(cp)]%*%cp
  bic.pred <- model.matrix(Salary~.,test)[,names(bic)]%*%bic
  cv.pred <- model.matrix(Salary~.,test)[,names(cv.coefs)]%*%cv.coefs
  # test errors 
  errors <- mean((r2.pred - test$Salary)**2)
  errors <- c(errors,mean((adjr2.pred - test$Salary)**2))
  errors <- c(errors,mean((cp.pred - test$Salary)**2))
  errors <- c(errors,mean((bic.pred - test$Salary)**2))
  errors <- c(errors,mean((cv.pred - test$Salary)**2))
  return(errors)
}


get.cv.error <- function(ncv, nmodels, data){
  # evaluates the mean cross-validation error of the linear model 
  # with the selected coefficients
  n.cv <- dim(data)[1]
  folds.cv <- cut(1:n.cv, ncv, labels=F)
  cv.errors <- matrix(nrow = ncv, ncol = nmodels)
  indices.cv <- 1:n.cv
  for(j in 1:ncv){
    test.indices.cv <- indices.cv[folds.cv==j]
    test.cv <- data[test.indices.cv,]
    train.cv <- data[-test.indices.cv,]
    cv.all.best<- regsubsets(x=Salary~.,data=train.cv,
                               nbest=1,nvmax=nmodels, # using all variables
                               method="forward" )
     for(m in 1:nmodels){
       cv.coefs <- coef(cv.all.best, id=m)
       cv.preds <- model.matrix(Salary~.,test)[,names(cv.coefs)]%*%cv.coefs
       # test errors 
       cv.errors[j,m] <- mean((cv.preds - test$Salary)**2)
    }
  }
  # selecting the model with the least mean error
  # expected test MSE estimated by CV for each model
  return(which.min(colMeans(cv.errors))) 
}

test.errors <- matrix(nrow=nfolds, ncol=5)

for(i in 1:nfolds){
  test.indices <- indices[folds==i]
  test <- Hitters[test.indices,]
  train <- Hitters[-test.indices,]
  # Now we'll use BSS on the train dataset 
  # And we'll record the error on the test set 
  # get best cv model 
  cv.best <- get.cv.error(ncv=5, nmodels=(dim(Hitters)[2]-1),data = train)
  test.errors[i,] <- get.bss.test.error(train=train, test=test, cv.best=cv.best)
  
}
```
Let's look at the results.
```{r}
test.errors <- data.frame(test.errors)
names(test.errors) <- c("r2","adjr2","cp","bic","cv")
test.errors
```
```{r}
plot(1:10, test.errors$r2, type="l", lty="dashed", col=2, ylab="test error", main="cv MSE estimate ", lwd=2)
lines(1:10, test.errors$adjr2, type="l", lty="dashed", col=3, lwd=2)
lines(1:10, test.errors$cp, type="l", lty="dashed", col=4, lwd=2)
lines(1:10, test.errors$bic, type="l", lty="dashed", col=5, lwd=2)
lines(1:10, test.errors$cv, type="l", lty="dashed", col=6, lwd=2)
legend("topright", legend = c("r2", "adjr2","cp","bic","cv"), col=c(2,3,4,5,6), lty="dashed")
```
```{r}
colMeans(test.errors)
which.min(colMeans(test.errors))
```
So the cross validation criteria seems to be the most reliable in model selection.
We'll now use this criteria to select the best model fitting it on the whole data. 
```{r}
best.cv <- get.cv.error(ncv=10, nmodels=(dim(Hitters)[2]-1), data=Hitters)
best.cv
```
Let's now look at the best model: 
```{r}
nmodels <- 19
all.best <- regsubsets(x=Salary~.,data=Hitters,
                       nbest=1,nvmax=nmodels, # using all variables
                       method="forward")
coef(all.best, id=9)
```
Let's look at some plots to see how the best model is seen according to other criteria. 
```{r}
par(mfrow=c(2,2))
s <- summary(all.best)
# rss
plot(s$rss,xlab="Model size",ylab="RSS",type="l")
points(9, s$rss[9], col="red",cex=2,pch=20)
# adjr2
plot(s$adjr2,xlab="Model size",ylab="Adjusted RSq",type="l")
points(9, s$adjr2[9], col="red",cex=2,pch=20)
# Cp
plot(s$cp,xlab="Model size",ylab="Cp",type='l')
points(9, s$cp[9], col="red", cex=2, pch=20)
# BIC
plot(s$bic,xlab="Model size",ylab="BIC",type='l')
points(9, s$bic[9], col="red", cex=2, pch=20)
```

```{r}
par(mfrow=c(1,1))
plot(all.best, scale="r2")
plot(all.best, scale="adjr2")
plot(all.best, scale="Cp")
plot(all.best, scale="bic")
```

## Shrinkage methods

We'll now loook at a different set of selection tools: shrinkage methods like LASSO and RIDGE. We'll actually use the ElasticNet model, of which Lasso and RIdge are special cases. 
```{r}
require(hdi)
data("riboflavin")
```
The riboflavin dataset records the riboflavin production by Bacillus subtilis together with the gene expressions. Each row refers to one gene, storing in x its expression level.
```{r}
require(glmnet)
attach(riboflavin)
lambda.grid <- grid <- 10^seq(10,-2, length = 100)
lasso <- glmnet(x = x, y=y, alpha = 1, lambda = lambda.grid)
ridge <- glmnet(x = x, y=y, alpha = 0, lambda = lambda.grid)
```
Exploration of the output: 
```{r}
dim(coef(lasso))
dim(coef(ridge))
```
```{r}
lasso$lambda[50]
round(coef(lasso)[,50],2)[coef(lasso)[,50]!=0]
ridge$lambda[50]
round(coef(ridge)[,50],2)[round(coef(ridge)[,50],2)!=0]
```
The above lambda is too high to allow any value to be different from zero. Let's have a broader look at the results with a plot: 
```{r}
plot(lasso)
```
```{r}
plot(ridge)
```
The above plots are a perfect synthesis of the differences between ridge and lasso: while ridge performs a "soft thresholding", slowly shrinking all variables to 0 as lambda increases, lasso performs a "hard thresholding", cutting off variables as they reach a certain threshold (determined by lambda).
```{r}
plot(ridge, xvar="lambda")
```
How to use a model for predictions? 
```{r}
preds <- predict(lasso, s = 0.01, newx=riboflavin$x)
mse.train <- mean((riboflavin$y - preds)**2)
mse.train

```


```{r}
preds.lasso <- predict(lasso, s = grid, newx=riboflavin$x)
plot(log(grid), colMeans((riboflavin$y - preds.lasso)**2), col="red", main="Train MSE", type="l", ylab = "", lty="dashed")

preds.ridge <- predict(ridge, s = grid, newx=riboflavin$x)
lines(log(grid), colMeans((riboflavin$y - preds.ridge)**2), col="blue", type="l", ylab = "", lty="dashed")

legend("right", legend=c("lasso","ridge"), col=c("red","blue"), lty="dashed")

```
Now let's use cross validation to do both model assessment and model selection: we're going to select the best lasso and ridge models on a train dataset and evaluate them with a second cross-validation against hold-out sets. 
```{r}
nfolds <- 10 
n <- dim(riboflavin)[1]
folds <- cut(1:n, breaks = nfolds, labels = F)
#shuffling
indices <- sample(1:n, size=n, replace = F)
res.cv <- matrix(nrow=10, ncol=4)
for(i in 1:nfolds){
  test.indices <- indices[folds==i]
  test <- riboflavin[test.indices, ]
  train <- riboflavin[-test.indices, ]
  ## model selection using cv on the train dataset
  ## LASSO
  lasso.cv <- cv.glmnet(x=train$x, y=train$y, 
                        alpha = 1, lambda = lambda.grid, nfolds =10)
  lasso.lambda <- lasso.cv$lambda.min
  res.cv[i,1] <- lasso.lambda
  lasso.fit <- glmnet(train$x, train$y, alpha =1, lambda = lasso.lambda)
  ## RIDGE
  ridge.cv <- cv.glmnet(x=train$x, y=train$y, 
                        alpha = 0, lambda = lambda.grid, nfolds = 10)
  ridge.lambda <- ridge.cv$lambda.min
  res.cv[i,3] <- ridge.lambda
  ridge.fit <- glmnet(train$x, train$y, alpha =0, lambda = ridge.lambda)
  ## model assessment on the test dataset
  ## LASSO
  lasso.predict <- predict(lasso.fit, newx=test$x)
  lasso.error<- mean((test$y-lasso.predict)**2)
  res.cv[i,2]<-lasso.error
  ## RIDGE
  ridge.predict <- predict(ridge.fit, newx=test$x)
  ridge.error<- mean((test$y-ridge.predict)**2)
  res.cv[i,4]<-ridge.error
}
res.cv <- data.frame(res.cv)
names(res.cv) <- c("lambda lasso","error lasso","lambda ridge", "error ridge")
```
```{r}
res.cv
```
```{r}
par(mfrow=c(2,2))
hist(res.cv$`error lasso`, main="LASSO error")
hist(res.cv$`lambda lasso`, main="LASSO lambda")
hist(res.cv$`error ridge`, main="RIDGE error")
hist(res.cv$`lambda ridge`, main="RIDGE lamdba")
```
Not only the ridge presents a higher expected error on the test set, but it also has more variability in the lambda values. Therefore we're going to use Lasso, fitting it to the whole dataset.
```{r}
lasso.cv <- cv.glmnet(x=x, y=y,alpha = 1, lambda = lambda.grid, nfolds =10)
plot(lasso.cv)
```
```{r}
best.lambda <- lasso.cv$lambda.min
lasso.fit <- glmnet(x,y, alpha=1, lambda=best.lambda, thresh=1e-12)
sum(coef(lasso.fit)!=0)
```
In the final model only 41 of the initial >4000 genes are active. 

# Selection bias

We'll now investigate the most common (and overlooked) mistake that comes with model selection.

Let's start by generating some data. 
```{r}
n <- 200
sigma <- 1
intercept <- 2
y <- intercept + rnorm(n=n, mean=0, sd=sigma)
X <- matrix(rnorm(n*7, mean=0, sd=1), nrow=n, ncol=7)
data <- data.frame(Y=y, X=X)
head(data)
pairs(data)
```
So, there's no relationship between Y and X, but let's suppose we don't know it and to make it more realistic let's give some fancy names to the X variables. 
```{r}
names(data)<-c("HealthIdx","Poverty12mo","MedianIncome","`%Obese"," InjuryRate","HeartRisk", "noConvict", "FamilyIssue")
head(data)
```
And now let's use our super-powerful model selection tools to look for the best way to model this data. 
```{r}
library(leaps)
p <- dim(data)[2]
bss.res <- regsubsets(HealthIdx~., data=data, nbest = 1, nvmax =p)
s <- summary(bss.res)
coef(bss.res, id=which.min(s$bic))
```
Okay our Best subset selection is telling us noCOnvict is the best predictor for our health index. Let's look at its p-value by fitting a linear model to it.
```{r}
fit.bss <- lm(HealthIdx~noConvict, data=data)
summary(fit.bss)
```
The p-value of no convict is not low enough, so no discovery! But will the p-value  always save us from false discoveries or was this just luck? 
But what if we repeat this approach multiple times? 
```{r}
nsim <- 1000
# we're going to record all p-values of the selected variables for each simulation
p.values <- matrix(nrow=nsim, ncol=p)
p.values <- data.frame(p.values)
names(p.values)<- names(data)[-1]
create.formula<- function(selected){
  formula <- selected[1]
  if(length(selected)>1){
    for(i in 2:length(selected)){
    formula <- paste(formula, selected[i], sep="+")
    }
  }
  return(formula)
}

for(j in 1:nsim){
  ##Let's draw again some data from the noise
  data$HealthIdx <- intercept + rnorm(n=n, mean=0, sd=sigma)
  ##Model selection 
  bss.res <- regsubsets(HealthIdx~., data=data, nbest = 1, nvmax =p)
  s <- summary(bss.res)
  selected <- names(coef(bss.res, id=which.min(s$cp)))[-1]
  ##Model assessment
  fit <- lm(paste("HealthIdx~",create.formula(selected)), data=data)
  s <- summary(fit)
  # extracting pvalues
  pvalues <- s$coefficients[selected,4]
  p.values[j, selected]=pvalues
}
```
```{r}
head(p.values)
```
Let's have a look at our results!
```{r}
hist(p.values$HeartRisk, main="HeartRIsk p-values after selection", breaks=20)
```
Let's now look at more honest p-values: the ones we get by fitting on the whole model, with no model selection.
```{r}
nsim <- 1000
# we're going to record all p-values of the selected variables for each simulation
p.values.noselection <- matrix(nrow=nsim, ncol=p)
p.values.noselection <- data.frame(p.values.noselection)
names(p.values.noselection)<- names(data)[-1]

for(j in 1:nsim){
  ##Let's draw again some data from the noise
  data$HealthIdx <- intercept + rnorm(n=n, mean=0, sd=sigma)
  # No model selection this time
  selected <- names(data)[-1] # all but the response
  fit <- lm("HealthIdx~.",data=data)
  s <- summary(fit)
  # extracting pvalues
  pvalues <- s$coefficients[-1,4]
  p.values.noselection[j, selected]=pvalues
}

```

```{r}
hist(p.values.noselection$HeartRisk, main="HeartRIsk p-values with no selection", breaks=20)

```
Now this is the true distribution of the p-values under the null. Why didn't we get this same distribution above? Because we applied what is called the *selection bias* to our analysis. Any time we use the data to make a decision (e.g. pick one model
instead of some others), we introduce a selection effect (bias). What is wrong with the naive approach we've seen first is that it's not taking into account that the test we're conducting is *conditioned on* the fact that that specific model has already been selected. 

Forward stepwise, Lasso, elastic net with cross-validation, etc, all use
the data in a way that would result in such bias.
Significance tests, prediction error, R2, goodness of fit tests, etc, will all suffer from selection bias.

So, how do we solve this?

### Solutions

The idea is basically to account for the conditioning, or, put in another way: if a variable “surprises” us enough to be included in the model, it must
surprise us again in order to be declared significant.

The first and easiest solution is to do what we've already done above: fit the whole model and look at those p-values. However, we should take into consideration a multiple testing issue and adjust our tests accordingly, since we're testing for the 0.5% on all the predictors. 

Another solution is to simply split the data in train and test set and NEVER use the test set before we've completed all the tests, selections and fitting. 
Let's look at the results of this approach on the above experiment: 
```{r}
nsim <- 1000
# we're going to record all p-values of the selected variables for each simulation
p.values.split <- matrix(nrow=nsim, ncol=p)
p.values.split <- data.frame(p.values.split)
names(p.values.split)<- names(data)[-1]

#we'll remove the noise that comes from the splitting and pre-determine the splitting before the simulation
train <- sample(c(TRUE,FALSE), size=(n*3/4), replace = T)
test <- (!train)

for(j in 1:nsim){
  ##Let's draw again some data from the noise
  data$HealthIdx <- intercept + rnorm(n=n, mean=0, sd=sigma)
  ##Model selection on the training set
  bss.res <- regsubsets(HealthIdx~., data=data[train,], nbest = 1, nvmax=p)
  s <- summary(bss.res)
  selected <- names(coef(bss.res, id=which.min(s$cp)))[-1]
  ##Model assessment on the test set
  fit <- lm(paste("HealthIdx~",create.formula(selected)), data=data[test,])
  s <- summary(fit)
  # extracting pvalues
  pvalues <- s$coefficients[selected,4]
  p.values.split[j, selected]=pvalues
}
```
```{r}
hist(p.values.split$HeartRisk, main="HeartRIsk p-values with selection and splitting", breaks=20)
```
Again we obtain honest p-values, which confirms the train-test splitting as a valid and selection-bias-free procedure. However, for how simple and robust (assumptions-wise) the splitting approach might be it has some drawbacks, the main one is the lack of reproducibility, due to the randomness introduced by the splitting. 

Some new research is working on *selective error control*, and here are some useful slides if you want to know more about it: http://joshualoftus.com/turing/shorttalk.pdf .