a) Training and test MSE [1 point]

Figure 2

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?

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?

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

b) Bias-variance trade-off [1 point]

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)?

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

Figure 5

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.

a) Understanding model output [1 point]

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

b) Model fit [1 point]

• 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?

# 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))

c) Confidence interval and hypothesis test [1 points]

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.

d) Prediction [1 point]

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.

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.

(Comment: Is that allowed - to only do the inverse? Yes, that could be the result of a first order Taylor expansion approximation. If you think this is interesting you can read more from page 36 here: https://www.math.ntnu.no/emner/TMA4267/2017v/L12.pdf - but that is not on the reading list of this course. You have probably used this result in your physics courses.)

• 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.

a) Logistic regression [1 point]

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\). (Hint: odds.)
• 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? (Hint: write out the equation for the boundary as \(x_2=intercept+slope*x_1\))
• 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.

Hint: if you have your calulated the slope and the intercept, your plot of train and test with boundary can be produced with:

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)+ggtitle("train and test and logistic boundary")

• 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?

Hint: if you have the predicted probabilities for class 1 for the test set in pred then you may do this to produce the confusion matrix:

testclass=ifelse(pred > 0.5, 1, 0)
t = table(test$y, testclass)
t

b) K-nearest neighbor classifier [0.5 point]

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?

Hint: you may not need this before d) if you just use the classification (and not the probabilites) when using KNN on the test set, but there is a slight complication on the output from prediction with KNN, since the probability of the winning class is given. Here is some code to avoid complications.

KNN3 = knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = F)
KNN3probwinning = attributes(knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = TRUE))$prob
KNN3prob <- ifelse(KNN3 == "0", 1-KNN3probwinning, KNN3probwinning)
# now KNN3prob is the probability of class 1 - and can be used for roc in d)

c) LDA (& QDA) [1.5 point]

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}\]
• We use the LDA 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.
• Perform LDA on the training data (using R). Just show the R-code for doing this
• Make the confusion table for the test set when using 0.5 as cut-off (Hint: just use predict on your lda-object). 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?

d) Compare classifiers [1 point]

• 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).

Hint: here is the code - you focus on interpretation!

library(pROC)
glmroc=roc(response=test$y,predictor=predglm)
plot(glmroc)
auc(glmroc)
KNN3 = knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = F)
KNN3probwinning = attributes(knn(train = train[,-1], test = test[,-1], k = 3, cl = train$y, prob = TRUE))$prob
KNN3prob <- ifelse(KNN3 == "0", 1-KNN3probwinning, KNN3probwinning)
KNN3roc=roc(response=test$y,predictor=KNN3prob) 
plot(KNN3roc)
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)
auc(lroc)