---
subtitle: "TMA4268 Statistical Learning V2018"
title: "Short solutions to Compulsory exercise 2"
author: "Thea Roksvåg, thea.roksvag@ntnu.no and Mette Langaas"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: #3rd letter intentation hierarchy
  # html_document:
  #   toc: true
  #   toc_float: true
  #   fig_caption: true
  pdf_document:
  #editor_options: 
  #  chunk_output_type: console
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,results="hold",message = FALSE,warning=FALSE)
````

Last changes: (23.05: corrected Q13 on acceleration/cylinder from Bjørn, 20.05 corrected Q2 decrease-increase from Rasmus.)

# Problem 1: Model selection and cross-validation

##1a) Theory
  
**Q1: How many possible models?** 
$2^d$ possible models: each covariate can either be in or out (2 possibilities) and we look at $d$ covariates.

**Q2: How to choose best model**

Use MSE to find the best model for each model complexity. Then choose the best model complexity from minimum BIC. 

If $R^2$ is used the most complex model will always be selected. Remember that RSS on the training data cannot increase when new covariates are added to the model, even if the covariates are random noise. This means that $R^2=1-\text{RSS}/\text{TSS}$ cannot decrease.

```{r,echo=FALSE,results="hide"}
library(ISLR)
ourAuto=data.frame("mpg"=Auto$mpg,"cylinders"=factor(cut(Auto$cylinders,2)),
                   "displace"=Auto$displacement,"horsepower"=Auto$horsepower,
                   "weight"=Auto$weight,"acceleration"=Auto$acceleration, 
                   "year"=Auto$year,"origin"=as.factor(Auto$origin))
colnames(ourAuto)
ntot=dim(ourAuto)[1]
ntot
set.seed(4268)
testids=sort(sample(1:ntot,ceiling(0.2*ntot),replace=FALSE))
ourAutoTrain=ourAuto[-testids,]
ourAutoTest=ourAuto[testids,]
library(leaps)
res=regsubsets(mpg~.,nbest=1,data=ourAutoTrain)
sumres=summary(res)
sumres
plot(res,scale="bic")
sumres$bic
which.min(sumres$bic)
coef(res,id=which.min(sumres$bic))
```

##1b) Interpreting output

**Q3: How to compare models**

To compare models with the same model complexity (same number of covariate) we may use the RSS or the $R^2$. The best model with two covariates is `mpg~weight+year`.

```{r, eval=TRUE,echo=TRUE}
fit=lm(mpg~.-acceleration,data=ourAutoTrain)
summary(fit)
library(nortest)
ad.test(rstudent(fit))
trainMSE=mean((fit$fitted.values-ourAutoTrain$mpg)^2)
trainMSE
pred=predict(fit,newdata=ourAutoTest)
Q5testMSE=mean((pred-ourAutoTest$mpg)^2)
Q5testMSE
```

**Q4: Choose model, BIC, best model:**

To compare models with different model complexity we could use penalized versions of the RSS, for example AIC, BIC and Mallow's Cp, or $R^2$adjusted - or use cross-validation to report MSE on test data. The smallest BIC is for the model with 7 coded covariates, `mpg~cylinders(5.5,8.01]+ displace+horsepower+weight+year+origin2+origin3`. This models explains 83% of the variability in the data, and all covariates are significant at level 0.004. However, normality of the residuals is questionable. MSE on training set is `r trainMSE`.

**Q5: Use this model fit to predict new values for `ourAutoTest` and report the MSE.**

 MSE on test set is `r Q5testMSE`.

##1c) Cross-validation

**Q6. Explain how $k$-fold cross-validation is performed. **

Divide training part of data into $k$ random divisions. Repeat for $j=1,...,k$: fit model(s) on all parts except $j$ and then use part $j$ for prediction - and then calculate a loss function. Sum over all parts and get a final MSE on the data left out. Use this to choose between different models - or to evaluate MSE.

For $k=5$ this is in detail: we divide the training data randomly into 5 folds of size $n/5$ each and call the folds $j=1$, to $j=5$. For $j=1,\ldots,5$:

* use the $4n/5$ observations from the folds except fold $j$ to fit the model
* the observations in the $j$th fold is left out and is the validation set, there are $n/5$ observations - and we denote them $(x_{0j},y_{0j})$, 
* we then predict $\hat{f}(x_{0j})$. 
* We calculate the error in the $j$th fold of the validation set as $\sum_{j} (y_{0j}-\hat{f}(x_{0j}))^2$ where the $j$ is for the validation fold

The total error on the validation set is thus the validation MSE=$\frac{1}{n}\sum_{m=1}^5 \sum_{j}(y_{0j}-\hat{f}(x_{0j}))^2$.

**Q7. Why may $k$-fold cross-validation be preferred to leave-one-out cross-validation?**

The LOOCV would give a less biased estimate (than $k$-fold) of the MSE on a test set, since the sample size used $(n-1)$ is close to the sample size to be used in the real world situation $(n)$, but the LOOCV is time consuming (however, for linear regression simple formula exists) and will be variable (the variance of the validation MSE is high). On the other hand the k-fold CV would give a biased estimate of the MSE on a test set since the sample size used is (k-1)/k of the sample size that would be used in the real world situation. However, the variance will be low and the procedure will be less time consuming. 


##1d) Programming cross-validation

**Q8. R-code for $10$-fold CV.**

```{r,eval=TRUE,echo=TRUE,results="hold"}
library(caret)
library(leaps)

nmodels=8
nfolds=10
set.seed(4268)
folds <- createFolds(ourAutoTrain$mpg,k=nfolds)
folds[[1]] #just checking which observations is in fold 1
folds[[2]]
lapply(folds,length)
sum(unlist(lapply(folds,length)))

# smart to make predict.regsubset
predict.regsubsets=function(object,newdata,id,...){
form=as.formula(object$call[[2]])
mat=model.matrix(form,newdata)
coefi=coef(object,id=id)
xvars=names(coefi)
mat[,xvars]%*%coefi
}

cv.errors=rep(0,nmodels)
# cross validation loop
for(j in 1:nfolds)
{
  thisfit=regsubsets(mpg~.,data=ourAutoTrain[-folds[[j]],],nvmax=nmodels)
  for(i in 1:nmodels)
  {
  pred=predict(thisfit,ourAutoTrain[folds[[j]],],id=i)
  cv.errors[i]=cv.errors[i]+sum((ourAutoTrain$mpg[folds[[j]]]-pred)^2)
  }
}
which.min(cv.errors)
# all covariates in - is the best

cvbest=lm(mpg~.,data=ourAutoTrain)
summary(cvbest)
```


**Q9. What is the optimal model complexity (number of parameters) in your regression?**

The optimal choice is $7$ parameters. NB: this is with seed 4268.

```{r}
# MSE on test set
pred=predict(cvbest,newdata=ourAutoTest)
Q10testMSE=mean((pred-ourAutoTest$mpg)^2)
Q10testMSE
Q5testMSE
```

**Q10. Report on the best model.**

Same model as in Q4 was chosen, so see Q5 for results.

# Problem 2: Shrinkage methods

##2a) Lasso and ridge regression

**Q11: Which figure is lasso and which is ridge.**

The $L_1$ penalty forces some of the coefficients to be zero. This is the case for the predictors in Figure 1, so this figure corresponds to the lasso. The $L_2$ penalty shrinks the coefficients to zero more smoothly. This is the case for the predictors in Figure 2, so this figure corresponds to ridge. 

**Q12: Explain effect of tuning parameter.**

The tuning parameter puts a penalty on the coefficient estimates. As the tuning parameter increases, the flexibility of the regression fit decreases because the coefficients are shrinked towards zero. The variance decreases, but the bias increases. $\lambda=0:$ Least squares solution.  $\lambda \rightarrow\infty:$ A model with intercept only. 

**Q13: Model selection?**

In ridge regression, all coefficients will shrink towards zero, but it will not set any of them exactly to zero (unless $\lambda=\infty$). Consequently, ridge regression will always generate a model involving all predictors. The lasso, however, forces some of the coefficient estimates to be exactly equal to zero when $\lambda$ is sufficiently large. Hence, the lasso is more suitable for performing variable selection and is more similar to best subset selection. 

In Figure 1 and Figure 2 we see that year and cylinders are the most important predictors. In 1b, we got that weight was the most important predictor, but year and cylinders were included in the best model with 3 covariates. 

```{r,echo=FALSE,results="hide"}
library(glmnet)
set.seed(4268)

x=model.matrix(mpg~.,ourAutoTrain)[,-1] #-1 to remove the intercept.
head(x)
y=ourAutoTrain$mpg

lambda=c(seq(from=5,to=0.1,length.out=150),0.01,0.0001) #Create a set of tuning parameters, adding low value to also see least squares fit
cv.out=cv.glmnet(x,y,alpha=1,nfolds=10,lambda=lambda, standardize=TRUE) #alpha=1 gives lasso, alpha=0 gives ridge

plot(cv.out)

```

##2b) Finding the optimal $\lambda$

**Q14: What does `cv.glmnet`?**

The function `cv.glmnet` uses cross validation to find the optimal $\lambda$, in this case 10-fold. Mean and standard deviation for validation set MSE is reported in a plot for each grid value of $\lambda$.

**Q15: What do you see in the plot?**

The plot shows the MSE with error estimates as a function of log(Lambda). The numbers on the top of the plot show the number of non-zero coefficients in the fitted model for each value of $\lambda$. One solution for optimal $\lambda$ is to choose the one where  MSE is at its minimum. 

Our textbook authors suggests a slight variation to this solution, called `lambda.1se`, find the $\lambda$ for the minimum MSE and then look at the standard deviation of the MSE at this point - and then move towards a more parsimonious model (increase $\lambda$) with mean MSE as the minimum MSE + one standard deviation. Read moe under `help(cv.glmnet)`.

**Q16: Finding the optimal lambda:**

```{r}
lambdaopt=cv.out$lambda.1se
print(lambdaopt)
```

##3c) Prediction

**Q17: Model fit**

Fit model, coefficients, ..
```{r}
lassofit=glmnet(x,y,alpha=1,lambda=lambdaopt,standardize=TRUE)
beta=coef(lassofit)
print(beta)
coef(lassofit)
```

**Q18: Prediction**

```{r}
# 0 for cylinder, displace, horsepower, weight, acceleration, year, 1 for origin2 and 0 for origin3
newx=matrix(c(0,150,100,3000,10,82,1,0),nrow=1)
pred_x0=predict(lassofit,s=lambdaopt,newx=newx)
# or - alternatively
sum(c(1,newx)*coef(lassofit))
```

The predicted `mpg` `r pred_x0`.

# Problem 3: Additiva non-linear regression

##3a) Additive model

**Q19: Fit additive model and comment.**

```{r}
library(gam)
fitgam=gam(mpg~bs(displace, knots=c(290))+
poly(horsepower,2)+weight+s(acceleration,df=3)+
origin,data=ourAutoTrain)
par(mfrow=c(2,3))
plot(fitgam,se=TRUE,col="blue")
summary(fitgam)
```

We see `displace` has two peaks, `horsepower` has the smallest CI for low values, the linear function in `weight` is very variable for small and high values, `acceleration` looks rather like a cubic function and there is a clear effect of `origin`.

**Q20: Basis for the cubic spline**

The cubic spline (see Module 7 page): Here $M=4$ for a cubic spline, and the basis is $x$, $x^2$,$x^3$, $(x-290)^3_{+}$. In `gam` an orthogonal version of this is used.

```{r}
M=4
cp=290
l=range(ourAutoTrain$displace)
par(mfrow=c(2,2),mar=c(4,1,2,1))
s=sapply(1:(M-1), function(y) curve(x^y, l[1],l[2],yaxt="n",
                                    xlab="displace",ylab="",main=bquote(b[.(y)](x)==x^.(y))))
s=sapply(1:1, function(y) curve(pmax(0,x-cp)^3,l[1],l[2],yaxt="n",
                                    xlab="displace",ylab="",main=bquote(b[.(y+M-1)](x-.(cp[y])["+"]^.(M-1)))))
```