---
subtitle: "TMA4268 Statistical Learning V2018"
title: "Compulsory exercise 1 with short solutions"
author: "Martina Hall, Martina.Hall@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
  pdf_document:
#editor_options: 
#  chunk_output_type: console
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(results = "hold")
knitr::opts_chunk$set(error=FALSE)
knitr::opts_chunk$set(warnings=FALSE)
```


# Problem 1 - Core concepts in statistical learning 

We consider a regression problem, where the true underlying curve is $f(x)=-x+x^2+x^3$ and we are considering $x \in [-3,3]$. 

This non-linear curve is only observed with added noise (either a random phenomenon, or unobservable variables influence the observations), that is, we observe $y=f(x)+\varepsilon$. In our example the error is sampled from $\varepsilon\sim N(0,2^2)$.

![Figure 1](Prob1f1.png)

In real life we are presented with a data set of pairs $(x_i,y_i)$, $i=1,\ldots,n$, and asked to provide a prediction at a value $x$. We will use the method of K nearest neighbour regression to do this here. 

We have a training set of $n=61$ observations $(x_i,y_i)$, $i=1,\ldots,n$. The KNN regression method provides a prediction at a value $x$ by finding the closes $K$ points and calculating the average of the observed $y$ values at these points. 

In addition we have a test set of $n=61$ observations (at the same grid points as for the training set), but now with new observed values $y$. 

We have considered $K=1,\ldots,25$ in the KNN method. Our experiment has been repeated $M=1000$ times (that is, $M$ versions of training and test set). 

## a) Training and test MSE 

![Figure 2](Prob1f2.png)

In the Figure 2 (above) you see the result of applying the KNN method with $K=1,2,10,25$ to our training data, repeated for $M$ different training sets (blue lines). The black lines show the true underlying curve. 

* Comment briefly on what you see.
* Does a high or low value of $K$ give the most flexible fit?

### Answers:
<span style="color:darkgreen">Here $K=1$ gives predicted values that on average (over the training sets) are not far from the true curve, but the large variability. Increasing the number of neighbours to $K=2$ gives less variance in the preditions, and still on average are not far from the true curve. Increasing the number of neighbours further to $K=10$, the variance is even more reduced, but the valued on the edge of our training set (at =-3 and x=3) are not well estimated, and $K=25$ gives little flexibility to the curve and doesn't fit well, especially on the edges of the data set. Hence, a low value of $K$ gives the most flexible fit.</span>


In Figure 3 (below) you see mean-squared errors (mean of squared differences between observed and fitted values) for the training set and for the test set (right panel for one training and one test set, and left panel for $M$). 

* Comment on what you see.
* What do you think is the "best" choice for K?


### Answers:
<span style="color:darkgreen">
We see that larger values of $K$ give larger variations in the MSE for both the test and the train set than for the smaller choices of $K$. For the training set, the MSE increases linearly with $K$, i.e. worse fit, but for the test set, there seems to be a small decrease at $K=5$ before it increases again. Based on this, we would choose $K=5$ as the best choice for $K$. </span>


Remark: in real life we do not know the true curve, and need to use the test data to decide on model flexibility (choosing $K$).

![Figure 3](Prob1f3.png)

## b) Bias-variance trade-off 

Now we leave the real world situation, and assume we know the truth (this is to focus on bias-variance trade-off). You will not observe these curves in real life - but the understanding of the bias-variance trade-off is a core skill in this course!

In the Figure 4 (below) you see a plot of estimated squared bias, estimated variance, true irreducible error and the sum of these (labelled total) and averaged over all values of $x$

The the squared bias and the variance is calculated based on the predicted values and the "true" values (without the added noise) at each $x$.  

* Explain how that is done. Hint: this is what the $M$ repeated training data sets are used for. 
* Focus on Figure 4. As the flexibility of the model increases ($K$ decreases), what happens with
    + the squared bias,  
    + the variance, and  
    + the irreducible error?
* What would you recommend is the optimal value of $K$? Is this in agreement with what you found in a)?



### Answers:


* <span style="color:darkgreen">How this is done?
The bias and variance of $\hat{f(x_0)}$ - when the true value is $f(x_0)$ is defined by
$$
\text{Bias}(\hat f(x_0)) =\text{E}[\hat f(x_0) - f(x_0)]
$$
$$
\text{Var}(\hat f(x_0)) = \text{E}[\hat f(x_0)^2] - \text{E}[\hat f(x_0)]^2
$$
Using the M resamples of the test data, we estimate the mean of  $\hat f(x_0)$ by the average over the $M$ predicted values at $x_0$, $\text{Ave}(\hat{f}(x_0))$, and the variance by the emipircal variance $\frac{1}{M-1} \sum_{i=1}^M (\hat{f}(x_0)-\text{Ave}(\hat{f}(x_0)))^2$.</span>
$$
\text{E}\left[\left(Y-\hat{f}\left(x_0\right)\right)^2\right]=\underbrace{\text{E}(\epsilon)}_{\text{Irredusible error}}+\underbrace{\text{Var}\left(\hat{f}\left(x_0\right)\right)}_{\text{Variance}}+\underbrace{\left[\text{E}\left(\hat{f}\left(x_0\right)\right)-\text{E}\left(Y\right)\right]^2}_{\text{Squared bias}}
$$

* <span style="color:darkgreen">
From Figure 4, we se that as K increases the squared bias increases rapidly, the variance is slightly reduced, but the irreducible error stays constant. Looking at the total MSE  (squared bias plus variance) we see that the optimal K-value lies somewhere between 3 and 5.</span>
* <span style="color:darkgreen">
Extra: We observe that the "optimal" $K$ varies for different $x_0$s. For $x_0=-2$ the best $K$ is between 5 and 10. For $x_0=0$ most values of $K$ is ok, and for $x_0=1$ the best is around $K=10$. For $x_0=2.5$ $K$ needs to be below 15. The reason for showing this is to see some of the $x_0$ values in our interval $[-3,3]$ that is given equal weight in deciding the optimal $K$.</span>



Extra: We have chosen to also plot curves at four values of $x$ - Figure 5 (below). Based on these four curves, that would you recommend is the optimal value of $K$? Is this in agreement with what you found previously (averaged over $x$)?

![Figure 4](Prob1f4.png)

![Figure 5](Prob1f5.png)

---

For completeness the R code used is given in the end of this file (listed here with `M=100` but `M=1000` was used). You do not need to run the code, this is just if you have questions about how this was done.


# Problem 2 - Linear regression 

The Framingham Heart Study is a study of the etiology (i.e. underlying causes) of cardiovascular
disease, with participants from the community of Framingham in Massachusetts, USA. 
For more more information about the Framingham Heart Study visit <https://www.framinghamheartstudy.org/>. The dataset used in here is subset of a teaching version of the Framingham data, used with permission from the Framingham Heart Study.

We will focus on modelling systolic blood pressure using data from n = 2600 
persons. For each person in the data set we have measurements of the
seven variables

* `SYSBP` systolic blood pressure,
* `SEX` 1=male, 2=female,
* `AGE` age in years at examination,
* `CURSMOKE` current cigarette smoking at examination: 0=not current smoker, 1= current smoker,
* `BMI` body mass index,
* `TOTCHOL` serum total cholesterol, and
* `BPMEDS` use of anti-hypertensive medication at examination: 0=not currently using, 1=currently using.

A multiple normal linear regression model was fitted to the data set with `-1/sqrt(SYSBP)` as
response and all the other variables as covariates. 

```{r,echo=TRUE,eval=TRUE}
library(ggplot2)
#data = read.table("https://www.math.ntnu.no/emner/TMA4268/2018v/data/SYSBPreg3uid.txt")
data = read.table("~/WWWemner/TMA4268/2018v/data/SYSBPreg3uid.txt")
dim(data)
colnames(data)
modelA=lm(-1/sqrt(SYSBP) ~ .,data = data)
summary(modelA)
```

## a) Understanding model output

We name the model fitted above `modelA`.

* Write down the equation for the fitted `modelA`.
* Explain (with words and formula) what the following in the `summary`-output means.
   + `Estimate` - in particular interpretation of `Intercept`
   + `Std.Error`
   + `t value`
   + `Pr(>|t|)`
   + `Residual standard error`
   + `F-statistic`



### Answers:


* <span style="color:darkgreen">Model A: 
$$-1/\sqrt{\text{SYSBP}}=\beta_0 + \beta_1 \text{SEX} + \beta_2 \text{AGE} + \beta_3 \text{CURSMOKE} + \beta_4 \text{BMI} + \beta_5 \text{TOTCHOL} + \beta_6 \text{BPMEDS} + \epsilon$$ with the fitted version
$$\widehat{1/\sqrt{\text{SYSBP}}}=-0.110 -0.0003 \text{SEX} +0.0002 \text{AGE} -0.0003 \text{CURSMOKE} + 0.0003 \text{BMI} + 0.00001 \text{TOTCHOL} + 0.0055 \text{BPMEDS}$$</span>

* <span style="color:darkgreen">The `Estimate` is the estimated regression coefficients, and are given by $\hat{\beta}=({\bf X}^T{\bf X})^{-1}{\bf X}^T{\bf Y}$. The interpretation of $\hat{\beta}_j$ is that when all other covariates are kept constant and the covariate $x_j$ is increased to from $x_j$ to $x_j+1$ then the response increases by $\hat{\beta}_j$. Example, holding all other variables constant, an increase of BMI from 25 to 26 will increase the response $-1/\sqrt{\text{SYSBP}}$ by $0.00031$. Similarily, for the binary variables, the coefficient estimates represents the change in the response when changing levels of the variable with one unit. For a female, the response will hence be reduced by $0.0003$ compared to a male (with the same values of all the other covariate). For all variables, negative value of the estimates give reduced response when increasing the corresponding variable, while positive estimates give increased response when increasing the corresponding variable.
The intercept, $\beta_0$ can be found by setting all other coefficients to zero. This involves also setting the covariate SEX to 0 - which has no meaning since SEX is coded as 1 for male and 2 for female.
</span>
* <span style="color:darkgreen">The `Std.Error` $\hat{SD}(\hat\beta_j)$ of the estimated coefficients is given by the square root of the diagonal entries of $({\bf X}^T{\bf X})^{-1}\hat{\sigma}^2$, where $\hat{\sigma}=\text{RSS}/(n-p-1)$. Here $n=2600$ and $p=6$.</span>
* <span style="color:darkgreen">The `t value` is the t-statistic $t = \frac{\hat\beta_j-\beta_j}{\hat{SD}(\hat\beta_j)}$, when we assume that $\beta_j=0$.</span>
* <span style="color:darkgreen">The `Pr(>|t|)` is the two-sided $p$-value for the null hypothesis $\beta_j=0$. The $p$-value is calculated as the probability of observing a test staistics equal to $|t|$ or larger in absolute value, assuming that the null hypothesis is true. A $p$-value less than 0.05 is considered statistically significant at a 5\% significance level.  </span>
* <span style="color:darkgreen">The residual standard error is the estimate of the standard deviation of $\epsilon$, and is given by RSS/$(n-p-1)$ where RSS=$\sum_{i=1}^n (y_i-\hat{y}_i)^2$. </span>
* <span style="color:darkgreen">The `F-statistic` is used test the hypothesis that all regression coefficients are zero, 
\begin{align*}
  H_0: & \beta_1 = \beta_2 = \cdots = \beta_p = 0 \quad \text{vs} \\
  H_1: &\text{at least one $\beta$ is $\neq 0$} \\ 
\end{align*}
and is computed by
\begin{equation*}
  F = \frac{(TSS-RSS)/p}{RSS/(n-p-1)}
\end{equation*}
where $TSS= \sum_{i=1}^n(y_i-\bar y)^2$, $RSS =  \sum_{i=1}^n(y_i-\hat y_i)^2$, $n$ is the number of observations and $p$ is the number of covariates (and $p+1$ the number of estimated regression parameters). If the $p$-value is less than 0.05, we reject the hypothesis that there are no coefficients with effect on the outcome in the model.  
</span>


## b) Model fit

* What is the proportion of variability explained by the fitted `modelA`? Comment.
* Use diagnostic plots of "fitted values vs. standardized residuals"" and "QQ-plot of standardized residuals" (see code below) to assess the model fit.
* Now fit a model, call this `modelB`, with `SYSBP` as response, and the same covariates as for `modelA`. Would you prefer to use `modelA` or `modelB` when the aim is to make inference about the systolic blood pressure?

```{r,eval=FALSE}
# residuls vs fitted
ggplot(modelA, aes(.fitted, .resid)) + geom_point(pch = 21) + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_smooth(se = FALSE, col = "red", size = 0.5, method = "loess") + 
  labs(x = "Fitted values", y = "Residuals", title = "Fitted values vs. residuals", subtitle = deparse(modelA$call))

# qq-plot of residuals
ggplot(modelA, aes(sample = .stdresid)) +
  stat_qq(pch = 19) + 
  geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals", title = "Normal Q-Q", subtitle = deparse(modelA$call))

# normality test
library(nortest) 
ad.test(rstudent(modelA))
```

### Answers:

* <span style="color:darkgreen">The $R^2$ statistic gives the proportion of variance explained by the model. In this model, the proportion of variability in $Y=-1/\sqrt{\text{SYSBP}}$ explained by the data $X$ is 0.2494. Since the range of $R^2$ is from 0 to 1, where for 1 all the variance in the response is explained by the regression model, we observe a fairly low number and we would have prefered it higher. However, these are medical data with low signal-to-noise ratio.</span>
* <span style="color:darkgreen">Looking at the diagnostic plots, the model fit looks good. The fitted values vs residuals plot is nice with semingly random spread and the QQ-plot looks nice as the plotted values follows the normal line. In addition, the Anderson-Darling normality test does not reject the hypothesis of normality.</span>
* <span style="color:darkgreen">For model B, we no longer model $-1/\sqrt{\text{SYSBP}}$, but rather SYSBP. This makes interpreation easier. However, looking at the diagnostic plots, we see that the QQ-plot looks suspicious at the tails, and the Anderson-Darling test rejects the null hypothesis of normal distribution.</span>

```{r,eval=T, include =T}

modelB = lm(SYSBP ~ .,data = data)
summary(modelB)

library(ggplot2)
# residuls vs fitted
ggplot(modelB, aes(.fitted, .resid)) + geom_point(pch = 21) + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_smooth(se = FALSE, col = "red", size = 0.5, method = "loess") + 
  labs(x = "Fitted values", y = "Residuals", title = "Fitted values vs. residuals", subtitle = deparse(modelB$call))

# qq-plot of residuals
ggplot(modelB, aes(sample = .stdresid)) +
  stat_qq(pch = 19) + 
  geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals", title = "Normal Q-Q", subtitle = deparse(modelB$call))

# normality test
library(nortest) 
ad.test(rstudent(modelB))
```

## c) Confidence interval and hypothesis test

We use `modelA` and focus on addressing the association between BMI and the response.

* What is the estimate $\hat{\beta}_{\text{BMI}}$ (numerically)?
* Explain how to interpret the estimated coefficient $\hat{\beta}_{\text{BMI}}$.
* Construct a 99% confidence interval for $\beta_{\text{BMI}}$ (write out the formula and calculate the interval numerically). Explain what this interval tells you.
* From this confidence interval, is it possible for you know anything about the value of the $p$-value for the test $H_0: \beta_{\text{BMI}}=0$ vs. $H_1:\beta_{\text{BMI}} \neq 0$? Explain.

### Answers:

* <span style="color:darkgreen"> $\hat\beta = (X^TX)^{-1}X^TY$. From the summary output we find that $\hat{\beta}_{\text{BMI}} = 0.0003$. This is the average increase in $-1/sqrt(SYSBP)$ for a unit increase in BMI. Hence, keeping all other covariates fixed - having a BMI of 24 instead of 23, the value of $-1/sqrt(SYSBP)$ will on average increase with 0.0003.</font>
* <span style="color:darkgreen">For linear regression where the distribution of the estimated coefficients are assumed to follow a t-distribution, we have that the $(1-\alpha)100$\%-confidence interval is given by
\begin{align*}
  \hat \beta \pm t_{\alpha/2, df} SD(\hat\beta)
\end{align*}
For $\hat\beta_{BMI}$ the 99\% confidence interval is hence given by
\begin{align*}
  [\hat \beta_{BMI} - t_{0.005, n-p-1} SD(\hat\beta_{BMI}), \hat \beta_{BMI} + t_{0.005, n-p-1} SD(\hat\beta_{BMI})] \\
\end{align*}
This means that before we have collected the data this interval has a 99% chance of covering the true value of $\beta_{BMI}$. After the interval is made - now this is [0.00023, 0.00038] the the true value is either within the interval or not. But, colleting new data and making 99% CIs, then 99% of these will on average cover the true $\beta_{BMI}$.</font>
* <span style="color:darkgreen">Since the interval does not cover 0, we know that the p-value is less than 0.01.</font>
```{r, include = T, eval = T}
n = dim(data)[1]
p = dim(data)[2]-1
betahat=modelA$coefficients[5]
sdbetahat=summary(modelA)$coeff[5,2]
UCI = betahat + qt(0.005, df = n-p-1, lower.tail = F)*sdbetahat
LCI = betahat - qt(0.005, df = n-p-1, lower.tail = F)*sdbetahat
c(LCI, UCI)
```

## d) Prediction 

Consider a 56 year old man who is smoking. He is 1.75 meters tall and his weight is 89 kilograms. His serum total cholesterol is 200 mg/dl and he is not using anti-hypertensive medication. 
```{r,eval=TRUE}
names(data)
new=data.frame(SEX=1,AGE=56,CURSMOKE=1,BMI=89/1.75^2,TOTCHOL=200,BPMEDS=0)
```

* What is your best guess for his `-1/sqrt(SYSBP)`? To get a best guess for his `SYSBP` you may take the inverse function of `-1/sqrt` (this would be a first order Taylor expansion).
* Construct a 90% prediction interval for his systolic blood pressure `SYSBP`. Comment. Hint: first contruct values on the scale of the response `-1/sqrt(SYSBP)` and then transform the upper and lower limits of the prediction interval.
* Do you find this prediction interval useful? Comment.

### Answers:

<span style="color:darkgreen"> Find best guess by prediction, and 90% prediction interval.</span>

```{r,eval=T}
pred = predict(modelA, newdata = new)
pred
f.inv = function(x) 1/x^2
sys = f.inv(pred)
#pred. interval
f.ci = predict(modelA, newdata = new, level = 0.9, interval = "prediction")
f.ci
sys.ci = f.inv(f.ci)
sys.ci
```

This prediction interval is very large and doesn't really tell us much. A person with our characteristics on average has a 90\% chance of having a systolic blood pressure between 108 and 168, and looking at the table given in <http://www.heart.org/HEARTORG/Conditions/HighBloodPressure/KnowYourNumbers/Understanding-Blood-Pressure-Readings_UCM_301764_Article.jsp#.WnLqWOYo_AI>, we sse that this interval covers almost all the levels from normal to high blood pressure. It seems our model is better for inference than prediction.</span>

# Problem 3 - Classification 

In this problem, we use a `wine` dataset of chemical measurement of two variables, `Color_intensity` and `Alcalinity_of_ash`, on 130 wines from two cultivars in a region in Italy. 

The data set is a subset of a data set from <https://archive.ics.uci.edu/ml/datasets/Wine>, see that page for information of the source of the data.

Below you find code to read the data, plot the data and to divide the data into a training set and a test set. To get your own unique division please change the seed (where it says `set.seed(4268)` you change 4268 to your favorite number).

```{r making data set, include=FALSE,eval=FALSE}
wine = read.csv(file = "https://www.math.ntnu.no/emner/TMA4268/2018v/data/wine.csv", dec = ".", head=F)
colnames(wine) = c("class", "Alcohol", "Malic_acid", "Ash", "Alcalinity_of_ash", "Magnesium", "Total_phenols", "Flavanoids", "Nonflavanoid_phenols", "Proanthocyanins", "Color_intensity", "Hue", "OD280/OD315_of_diluted wines", "Proline")
ourwine = wine[which(wine$class!=3),c(1,5,11)]
write.table(ourwine,"Comp1Wine.csv")
```

```{r,warning=FALSE,results="hold",message=FALSE,error=FALSE}
library(ggplot2)
library(GGally)
library(class)
library(MASS)
library(pROC)
#wine=read.csv("https://www.math.ntnu.no/emner/TMA4268/2018v/data/Comp1Wine.csv",sep=" ")
wine = read.csv("~/WWWemner/TMA4268/2018v/data/Comp1Wine.csv",sep=" ")
wine$class=as.factor(wine$class-1)
colnames(wine)=c("y","x1","x2")
ggpairs(wine, ggplot2::aes(color=y))

n=dim(wine)[1]
set.seed(4268) #to get the same order if you rerun - but you change this to your favorite number
ord = sample(1:n) #shuffle 
test = wine[ord[1:(n/2)],]
train = wine[ord[((n/2)+1):n],]
```

In our data the two classes are named `y` and coded $Y=0$ and $Y=1$, and we name $x_1$= `Alcalinity_of_ash`=`x1` and $x_2$=`Color_intensity`=`x2`.

## a) Logistic regression 
We assume a logistic regression model for observation $i$, $i=1,\ldots,n$:
$$
\text{Pr}(Y_i = 1| {\bf X}={\bf x}_i) = p_i = \frac{e^{\beta_0 + \beta_1x_{i1} + \beta_2 x_{i2}}}{ 1+ e^{\beta_0 + \beta_1x_{i1} + \beta_2x_{i2}}}
$$

* Use this expression to show that logit$(p_i)=\log(\frac{p_i}{1-p_i})$ is a linear function.
* Fit a logistic regression model on `y~x1+x2` to the training set.
* Give an interpretation of $\hat{\beta}_1$ and $\hat{\beta}_2$.
* We use the rule to classify to class 1 for an observation with covariates ${\bf x}$ if $\hat{\text{Pr}}(Y=1\mid {\bf x})>0.5$. Write down the formula for the class boundary between the classes. What type of boundary is this?
* Make a plot with the training observations and the class boundary. Hint: in `ggplot` points are added with `geom_point` and a line with `geom_abline(slope=b, intercept=a)` where $a$ and $b$ comes from your class boundary, and title with `ggtitle`.
* Use the `summary` output to manually derive the predicted probability $\hat{\text{Pr}}(Y=1\mid x_1=17, x_2=3)$. What is the interpretation of this value?
* Compute predicted probabilites for all observations in the test set.
* Make the confusion table for the test set when using 0.5 as cutoff for the probabilities. Calculate the sensitivity and specificity on the test set. How would you evaluate the performance of this classification?

### Answers:

* <span style="color:darkgreen"> 
\begin{align*}
 \log\big( \frac{p_i}{1-p_i} \big) &= \log \Big( \frac{\frac{e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}{1+e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}}{1- \frac{e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}{1+ e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}} \Big) = \log \big(\frac{\frac{e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}{1+ e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}}{\frac{1}{1+ e^{\beta_0 + \beta_1x_1+ \beta_2x_2}}} \big) \\
 &= \log(e^{\beta_0 + \beta_1x_1+ \beta_2x_2}) = \beta_0 + \beta_1x_1+ \beta_2x_2 \\
\end{align*}
</span>
```{r, include = TRUE, eval=TRUE}
fit = glm(y~x1+x2, data = train, family = "binomial")
summary(fit)
```
* <span style="color:darkgreen"> $\hat{\beta}_1$ is the estimated coefficient for `x1` with a value 0.1332. This means that the odds of begin a wine of class 2 is multiplied by $\exp(\hat\beta_1)=1.14$ if `x1` increases by one unit.
Further, $\hat{\beta}_2$ is the estimated coefficient for `x2` with value -2.34, which is interpreted as the odds of begin a wine of class 2 is multiplied by $\exp(\hat\beta_2)=0.097$ if `x2` increases by one unit.
</span>
* <span style="color:darkgreen"> The class boundary is where there is equal probability of both classes, $p_i = 0.5$. Using the logit functon, we find that $\frac{p_i}{1-p_i} = \exp(\beta_0 + \beta_1 x_1 + \beta_2x_2)$. With $p_i = 0.5$, the class boundary is given by $1 = \exp(\beta_0 + \beta_1 x_1 + \beta_2x_2)$ and hence $0 = \beta_0 + \beta_1 x_1 + \beta_2x_2$, and $x_2=-\frac{\beta_0}{\beta_2}-\frac{\beta_1}{\beta_2}x_1$. The class boundary is a straight line and can easily be added to the scatter plot below.</span>

```{r, include = T, eval=T}
betas = fit$coefficients
slope = -betas[2]/betas[3]
intercept = -betas[1]/betas[3]

g1 = ggplot(data=train,aes(x=x1, y=x2, colour=y)) + geom_point(pch = 1)
g1 + geom_point(data = test, pch = 3) + geom_abline(slope=slope,intercept=intercept)
```

```{r, include = T, eval=T}
x=c(1,17,3)
p = exp(fit$coefficients%*%x)/(1+exp(fit$coefficients%*%x))
p
```
For a wine with `x1` of 17 and  `x2` of 3, there is a 93% probability that the wine is of class 1.

```{r, include=TRUE, eval=TRUE}
pred = predict.glm(fit, newdata = test, type = "response")
predglm = predict.glm(fit, newdata = test, type = "response")
testclass=ifelse(predglm > 0.5, 1, 0)
t = table(test$y, testclass)
t

n = length(test$y)
error = (n-sum(diag(t)))/n
error
library(caret)
confusionMatrix(as.factor(testclass),reference=as.factor(test$y),positive="1")
```
* <span style="color:darkgreen"> 10 misclassifications equally divided over both wines, and a error rate pf 0.15. I would say that this classification works pretty well. The sensitivity is $30/35=0.86$ (true class 1), and the specificity is $25/30=0.83$ (true class 0). This is given that we define class 1 as the "positive class" - if we do differently the sensitivity and specificity will be swapped. For a disease situation it is "easy" to choose the disease as coded as the "postive", but for our wine example - we may use either 0 or 1 as positive. So, if  you got sensitivity and specificity swapped - that is ok.
</span>

## b) K-nearest neighbor classifier 
To decide the class of an new observation, the KNN classifier uses the nearest neighbours in the following way,
$$
P(Y=1|X=x_0) = \frac{1}{K}\sum_{i\in \mathcal{N_0}}I(y_i=j).
$$

* Explain this expression does, and what the different elements are.
* Use KNN with $K=3$ to classify the wines in the test set. 
* Make the confusion table for the test set when using 0.5 as cutoff. Calculate the sensitivity and specificity on the test set. How would you evaluate the performance of this classification? 
* Repeat with $K=9$. Which of these two choices of $K$ would you prefer and why? Why don't we just choose $K$ as high or as low as possible? 

### Answers:

* <span style="color:darkgreen"> Given an integer $K$ and a test observation $x_0$, the KNN classifier first identifies the $K$ points in the training data that are closest to $x_0$, represented by $\mathcal{N_0}$. It then estimates the conditional probability for class $j$ as the fraction of points in $\mathcal{N_0}$ whose response values equal $j$. </span>

```{r, include=TRUE, eval=TRUE}
KNN3 = knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = F)
t3 = table(test$y, KNN3)
t3
apply(t3,1,sum)
n = length(test$y)
error = (n-sum(diag(t3)))/n
error
KNN3probwinning = attributes(knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = TRUE))$prob
KNN3prob <- ifelse(KNN3 == "0", 1-KNN3probwinning, KNN3probwinning)
#cbind(KNN3prob,KNN3,KNN3probwinning) to check that this is correct
KNN3roc=roc(response=test$y,predictor=KNN3prob)
cbind(KNN3roc$threshold, KNN3roc$sens, KNN3roc$specificities)
```

* <span style="color:darkgreen"> KNN3:
This gives a rather good result with a misclassification rate of 0.12. Since we have only three neighbours to base the classification on the only possible values for the probability of class 2 is 0, 0.33, 0.67 and 1. For a cut-off of 0.5 the sensitivity is 32/35=0.91 and the specificity is 25/30=0.83.
</span>


```{r, include=TRUE, eval=TRUE}
KNN9= knn(train = train[,-1], test = test[,-1], k = 9, cl = train$y, prob = F)
t9 = table(test$y, KNN9)
t9
error = (n-sum(diag(t9)))/n
error
KNN9probwinning = attributes(knn(train = train[,-1], test = test[,-1], k = 9, cl = train$y, prob = TRUE))$prob
KNN9prob <- ifelse(KNN9 == "0", 1-KNN9probwinning, KNN9probwinning)
KNN9roc=roc(response=test$y,predictor=KNN9prob)
cbind(KNN9roc$threshold, KNN9roc$sens, KNN9roc$specificities)
table(KNN3,KNN9)
```

* <span style="color:darkgreen"> KNN9
For $K=9$ the misclassification rate when using 0.5 as cut-off is 0.15, so higher than for $K=3$, but similar to logistic regression. The sensitivity is 
The classifications made by $K=3$ and $K=9$ are not that different, 27+34=61 common classifications. But, since $K=3$ have the lowest misclassification rate we might prefer that, however for $K=3$ we may have strange boundary effects.</span>

* <span style="color:darkgreen"> If we choose $K=1$ that might lead to a too flexible class boundary, and with $K=n$ this might be too inflexible.
</span>


## c) LDA (& QDA) 
In linear discriminant analysis, with $K$ classes, we assign a class to a new observation based on the posterior probability 

$$\text{Pr}(Y=k | {\bf X}={\bf x}) = \frac{\pi_k f_k({\bf x})}{\sum_{l=1}^K \pi_l f_l({\bf x})},$$
where 
$$f_k({\bf x}) = \frac{1}{(2 \pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}}e^{-\frac{1}{2}({\bf x}-\boldsymbol{\mu_k})^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol{\mu_k})}.$$

* Explain what is $\pi_k$, $\boldsymbol{\mu}_k$, $\boldsymbol{\Sigma}$ and $f_k(x)$ in our `wine` problem.
* How can we estimate $\pi_k$, $\boldsymbol{\mu}_k$ and $\boldsymbol{\Sigma}$? Compute estimates for these quantities based on the training set.

In a two class problem ($K=2$) the decision boundary for LDA between class 0 and class 1 is where $x$ satisfies
$$
\text{Pr}(Y=0 | {\bf X}={\bf x}) = \text{Pr}(Y=1 | {\bf X}={\bf x}).
$$

* Show that we can express this as 
\begin{align}
\delta_0({\bf x}) &= \delta_1({\bf x}),
\end{align}
where
\begin{align}
\delta_k({\bf x}) &= {\bf x}^T \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_k - \frac{1}{2}\boldsymbol{\mu}_k^T \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_k + \log \pi_k; \quad k\in\{0,1\}.
\end{align}
* Perform LDA on the training data (using R).
* We use the rule to classify to class 1 for an observation with covariates ${\bf x}$ if $\hat{\text{Pr}}(Y=1\mid {\bf x})>0.5$. Write down the formula for the class boundary between the classes.
* Make a plot with the training observations and the class boundary. Add the test observations to the plot (different markings). Hint: in `ggplot` points are added with `geom_points` and a line with `geom_abline(slope=b, intercept=a)` where $a$ and $b$ comes from your class boundary.
* Make the confusion table for the test set when using 0.5 as cut-off. Calculate the sensitivity and specificity on the test set. How would you evaluate the performance of this classification? 
* If you where to perform QDA instead of LDA, what would be the most important difference between the QDA and LDA philosophy?

### Answers:

* <span style="color:darkgreen"> Here $\pi_k$ is the prior probability that a randomly chosen observation comes from the $k$th class.
We assume that the observations of class $k$ comes from a multivariate normal distribution $f_k(x)$, where $\boldsymbol{\mu}_k$ is the mean of the $k$th class and $\Sigma$ is the variance.</span>
* <span style="color:darkgreen">We estimate these in the following way:
\begin{align*}
  \hat{\pi}_k &= n_k/n \\
  \hat{\mu}_k &= \frac{1}{n_k}\sum_{i: y_i=k}x_i \\ 
\hat{\boldsymbol{\Sigma}}_k&=\frac{1}{n_k-1}\sum_{i:y_i=k} ({\bf X}_i-\hat{\boldsymbol{\mu}}_k ) ({\bf X}_i-\hat{\boldsymbol{\mu}}_k)^T\\
\hat{\boldsymbol{\Sigma}}&= \sum_{k=1}^K \frac{n_k - 1}{n - K} \cdot \hat{\boldsymbol{\Sigma}}_k.
\end{align*}

For the training dataset we have that

```{r}
n=dim(train)[1]
train0 = train[which(train$y==0),2:3]
train1 = train[which(train$y==1),2:3]
print("pi0 and pi1")
pi0=dim(train0)[1]/n; pi1=dim(train1)[1]/n
c(pi0,pi1)
print("mu")
mu0=apply(train0,2,mean)
mu1=apply(train1,2,mean)
mu0
mu1
print("Sigma")
Sigma=((dim(train0)[1]-1)*var(train0)+(dim(train1)[1]-1)*var(train1))/(n-2)
Sigma
```
</span>

* <span style="color:darkgreen"> Show $\delta$:
\begin{align*}
 \text{Pr}(Y=0 | X={\bf x}) = \text{Pr}(Y=1 | X={\bf x}) &= \frac{\pi_0f_0(x)}{\pi_0f_0(x)+\pi_1f_1(x)} = \frac{\pi_1f_0(1)}{\pi_0f_0(x)+\pi_1f_1(x)} \\
\pi_0 f_0({\bf x}) = \pi_1 f_1({\bf x}) = \pi_0 f_0({\bf x}) = 
\pi_0\frac{1}{(2 \pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}}
e^{\frac{1}{2}({\bf x}-\boldsymbol{\mu_0})^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol{\mu_0})} &=  \pi_1\frac{1}{(2 \pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}}e^{\frac{1}{2}({\bf x}-\boldsymbol{\mu_1})^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol{\mu_1})} \\
\log(\pi_0) -\frac{1}{2}({\bf x}-\boldsymbol{\mu}_0)^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol{\mu}_0) &= \log(\pi_1) -\frac{1}{2}({\bf x}-\boldsymbol{\mu}_1)^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol{\mu}_1) \\
\log(\pi_0)  -\frac{1}{2}x^T\boldsymbol{\Sigma}^{-1}{\bf x} + {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_0 - \frac{1}{2}\mu_0^T\boldsymbol{\Sigma}^{-1}\mu_0 &= \log(\pi_1)  -\frac{1}{2}{\bf x}^T\boldsymbol{\Sigma}^{-1}{\bf x} + {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_1 - \frac{1}{2}\mu_1^T\boldsymbol{\Sigma}^{-1}\mu_1 \\
\log(\pi_0) + {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_0 - \frac{1}{2}\mu_0^T\boldsymbol{\Sigma}^{-1}\mu_0 &= \log(\pi_1) + {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_1 - \frac{1}{2}\mu_1^T\boldsymbol{\Sigma}^{-1}\mu_1 = \delta_0({\bf x}) = \delta_1({\bf x}) \\
\end{align*}
</span>
* <span style="color:darkgreen"> LDA in R, and draw class boundary. Have added also logistic regression boundary.</span>
```{r}
ltrain=lda(y~x1+x2,data=train)
a=solve(Sigma)%*%(mu0-mu1)
c=-0.5*t(mu0)%*%solve(Sigma)%*%mu0+0.5*t(mu1)%*%solve(Sigma)%*%mu1+log(pi0)-log(pi1)
interceptL = -c/a[2]
slopeL = -a[1]/a[2]

g1 = ggplot(data=train,aes(x=x1, y=x2, colour=y)) + geom_point(pch = 1)
g1 + geom_point(data = test, pch = 3) + geom_abline(slope=slope,intercept=intercept,colour="green")+ geom_abline(slope=slopeL,intercept=interceptL,colour="darkblue")

```

* <span style="color:darkgreen"> Formula for LDA class boundary (0.5 cut-off):
$\delta_0({\bf x}) = \delta_1({\bf x})$, and thus $\delta_0({\bf x})- \delta_1({\bf x})=0$.
$$\log(\pi_0) + {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_0 - \frac{1}{2}\mu_0^T\boldsymbol{\Sigma}^{-1}\mu_0  -\log(\pi_1) - {\bf x}^T\boldsymbol{\Sigma}^{-1}\mu_1 + \frac{1}{2}\mu_1^T\boldsymbol{\Sigma}^{-1}\mu_1 =0$$
$${\bf x}^T\boldsymbol{\Sigma}^{-1}(\mu_0-\mu_1) - \frac{1}{2}\mu_0^T\boldsymbol{\Sigma}^{-1}\mu_0 + \frac{1}{2}\mu_1^T\boldsymbol{\Sigma}^{-1}\mu_1 + \log(\pi_0)-\log(\pi_1)=0$$
Letting all the terms except the first be noted $c$, then $c$ is found in the R printout above.
The equation is then of the form
$ax_1+bx_2+c=0$, so $x_2=-c/b -a/b x_1$, as calculated above. This gave
```{r}
interceptL
slopeL
```
</span>

* <span style="color:darkgreen"> Confusion table, with misclassification rate 5+6/65=0.17, sensitivity 30/35=0.86 and specificity 24/30=0.8.
</span>

* <span style="color:darkgreen">If we went from LDA to QDA we would allow the classes to have differen covariance matrices. This would give nonlinear (quadratic) class boundaries. It is not clear if that would be a good choice. </span>

```{r}
ltrain=lda(y~x1+x2,data=train)
lpred=predict(object = ltrain, newdata = test)$posterior[,2]
lroc=roc(response=test$y,lpred)
testclass=ifelse(lpred > 0.5, 1, 0)
t = table(test$y, testclass)
t

n = length(test$y)
error = (n-sum(diag(t)))/n
error

```

## d) Compare classifiers 
* Compare your results from the different classification methods (logistic regression, your preferred KNN, LDA) based on the 0.5 cut-off on posterior probability classification rule. Which method would you prefer?
* Explain what an ROC curve is and why that is useful. Would your preference (to which method is the best for our data) change if a different cut-off was chosen?  Answer this by producing ROC-curves for the three methods on the test set. Also calculate AUC. Hint: use function `res=roc(response=test$y,predictor)` in `library(pROC)` where the predictor is a vector with your predicted posterior probabilites for the test set, and then `plot(res)` and `auc(res)`.

```{r}
glmroc=roc(response=test$y,predictor=predglm)
plot(glmroc)# logistic solid line
print("logistic regression")
auc(glmroc)
KNN3roc=roc(response=test$y,predictor=KNN3prob) #see above for code
plot(KNN3roc,add=TRUE,lty=2) # KNN3 dashed
print("KNN3")
auc(KNN3roc)
ltrain=lda(y~x1+x2,data=train)
lpred=predict(object = ltrain, newdata = test)$posterior[,1]
lroc=roc(response=test$y,lpred)
plot(lroc,add=TRUE,lty=3) # LDA dotted
print("LDA")
auc(lroc)
```

<span style="color:darkgreen">
Best method wrt AUC: LDA, but not very different from logistic regression, while 3KNN is best wrt misclassification rate on test set. 
</span>

<span style="color:darkgreen">
In the ROC-curve we plot the sensitivity and 1- spesificity against each other for all possible thresholds of the probability for class 1. To construct the ROC-curve we would have to calculate the sensitivity and spesificity for different values of the cutoff $p(x)>cut$. Using a threshold of 0.5, you say that if a new person has a probability of 0.51 of having the disease, he is classified as diseased. Another person with a probability of 0.49 would then be classified as non-diseased. The ROC-curve and the area under the ROC-curve are useful tools as they consider all possible thresholds for the cutoff. 
</span>

<span style="color:darkgreen">
The AUC is the area under the ROC-curve and gives the overall performance of the test for all possible thresholds. An AUC value of 1 means a perfect fit for all possible thresholds, while a AUC of 0.5 corresponds to the classifyer that performs no better than chance. Hence, a classification method $p(x)$ giving 0.6 and another $q(x)$ giving 0.7, we would prefer $q(x)$ as it has the highest AUC value - in general - unless there is a specific reason for only wanting to consider one specific value for the cut.off.
</span>

---

# R-code for Problem 1

\footnotesize
```{r,eval=FALSE,echo=TRUE}
library(FNN)
library(ggplot2)
library(ggpubr)
library(reshape2)
maxK=25
M=100 # repeated samplings, x fixed  - examples were run with M=1000
x = seq(-3, 3, 0.1)
dfx=data.frame(x=x)
truefunc=function(x) return(-x+x^2+x^3)
true_y = truefunc(x)

set.seed(2) # to reproduce
error = matrix(rnorm(length(x)*M, mean=0, sd=2),nrow=M,byrow=TRUE)
testerror = matrix(rnorm(length(x)*M, mean=0, sd=2),nrow=M,byrow=TRUE)
ymat = matrix(rep(true_y,M),byrow=T,nrow=M) + error
testymat = matrix(rep(true_y,M),byrow=T,nrow=M) + testerror

ggplot(data=data.frame(x=x,y=ymat[1,]),aes(x,y))+geom_point(col="purple",size=2)+stat_function(fun=truefunc,lwd=1.1,colour="black")+ggtitle("Training data")

predarray=array(NA,dim=c(M,length(x),maxK))
for (i in 1:M)
{
  for (j in 1:maxK)
  {
    predarray[i,,j]=knn.reg(train=dfx,test=dfx,y=c(ymat[i,]),k=j)$pred
  }
}
# first - just plot the fitted values - and add the true curve in black
# M curves and choose k=1,2,10,30 in KNN

# rearranging to get data frame that is useful
thislwd=1.3
stackmat=NULL
for (i in 1:M) stackmat=rbind(stackmat,cbind(x,rep(i,length(x)),predarray[i,,]))
colnames(stackmat)=c("x","rep",paste("K",1:maxK,sep=""))
sdf=as.data.frame(stackmat) 
yrange=range(apply(sdf,2,range)[,3:(maxK+2)])
# making the four selected plots
p1=ggplot(data=sdf,aes(x=x,y=K1,group=rep,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
p1=p1+stat_function(fun=truefunc,lwd=thislwd,colour="black")+ggtitle("K1")
p2=ggplot(data=sdf,aes(x=x,y=K2,group=rep,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
p2=p2+stat_function(fun=truefunc,lwd=thislwd,colour="black")+ggtitle("K2")
p10=ggplot(data=sdf,aes(x=x,y=K10,group=rep,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
p10=p10+stat_function(fun=truefunc,lwd=thislwd,colour="black")+ggtitle("K10")
p25=ggplot(data=sdf,aes(x=x,y=K25,group=rep,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
p25=p25+stat_function(fun=truefunc,lwd=thislwd,colour="black")+ggtitle("K30")
ggarrange(p1,p2,p10,p25)

# calculating trainMSE and testMSE
trainMSE=matrix(ncol=maxK,nrow=M)
for (i in 1:M) trainMSE[i,]=apply((predarray[i,,]-ymat[i,])^2,2,mean)
testMSE=matrix(ncol=maxK,nrow=M)
for (i in 1:M) testMSE[i,]=apply((predarray[i,,]-testymat[i,])^2,2,mean)
#rearranging to get data frame that is useful
stackmat=NULL
for (i in 1:M) stackmat=rbind(stackmat,cbind(rep(i,maxK),1:maxK,trainMSE[i,],testMSE[i,]))
colnames(stackmat)=c("rep","K","trainMSE","testMSE")
sdf=as.data.frame(stackmat) 
yrange=range(sdf[,3:4])
# plotting traning and test MSE
p1=ggplot(data=sdf[1:maxK,],aes(x=K,y=trainMSE))+scale_y_continuous(limits=yrange)+geom_line()
pall= ggplot(data=sdf,aes(x=K,group=rep,y=trainMSE,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
testp1=ggplot(data=sdf[1:maxK,],aes(x=K,y=testMSE))+scale_y_continuous(limits=yrange)+geom_line()
testpall= ggplot(data=sdf,aes(x=K,group=rep,y=testMSE,colour=rep))+scale_y_continuous(limits=yrange)+geom_line()
ggarrange(p1,pall,testp1,testpall)


# calculating bias^2 and variance
meanmat=matrix(ncol=length(x),nrow=maxK)
varmat=matrix(ncol=length(x),nrow=maxK)
for (j in 1:maxK)
{
  meanmat[j,]=apply(predarray[,,j],2,mean) # we now take the mean over the M simulations - to mimic E and Var at each x value and each KNN model
  varmat[j,]=apply(predarray[,,j],2,var)
}
bias2mat=(meanmat-matrix(rep(true_y,maxK),byrow=TRUE,nrow=maxK))^2 #here the truth is finally used!

# preparing to plot
df=data.frame(rep(x,each=maxK),rep(1:maxK,length(x)),c(bias2mat),c(varmat),rep(4,prod(dim(varmat)))) #irr is just 4
colnames(df)=c("x","K","bias2","variance","irreducible error") #suitable for plotting
df$total=df$bias2+df$variance+df$`irreducible error`
hdf=melt(df,id=c("x","K"))
# averaged over all x - to compare to train and test MSE
hdfmean=
  hdf %>%
  group_by(K,variable) %>%
  summarise (mean_value=mean(value))
ggplot(data=hdfmean[hdfmean[,1]<31,],aes(x=K,y=mean_value,colour=variable))+geom_line()+ggtitle("averaged over all x")

# extra: what about different values of x?
hdfatxa=hdf[hdf$x==-2,]
hdfatxb=hdf[hdf$x==0,]
hdfatxc=hdf[hdf$x==1,]
hdfatxd=hdf[hdf$x==2.5,]
pa=ggplot(data=hdfatxa,aes(x=K,y=value,colour=variable))+geom_line()+ggtitle("x0=-2")
pb=ggplot(data=hdfatxb,aes(x=K,y=value,colour=variable))+geom_line()+ggtitle("x0=0")
pc=ggplot(data=hdfatxc,aes(x=K,y=value,colour=variable))+geom_line()+ggtitle("x0=1")
pd=ggplot(data=hdfatxd,aes(x=K,y=value,colour=variable))+geom_line()+ggtitle("x0=2.5")
ggarrange(pa,pb,pc,pd)
```

\normalsize





