---
title: "Hwk #2: Classification methods and Penalization"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(MASS)
```

For this homework we will use NHANES data that exists in a package for R.

NHANES consists of survey data collected by the US National Center for Health Statistics (NCHS) which has conducted a series of health and nutrition surveys since the early 1960's. Since 1999 approximately 5,000 individuals of all ages are interviewed in their homes every year and complete the health examination component of the survey. The health examination is conducted in a mobile examination center (MEC).

Note that there is the following warning on the NHANES website:
“For NHANES datasets, the use of sampling weights and sample design variables is recommended for all analyses because the sample design is a clustered design and incorporates differential probabilities of selection. If you fail to account for the sampling parameters, you may obtain biased estimates and overstate significance levels.”

For this homework, please ignore this warning and just apply our analyses to the data as if they were randomly sampled! We will be using the data called `NHANESraw`.

For questions that ask for your comments, it suffices to answer with one or two sentences in each case.

## Data Preparation

1. Install the package `NHANES` into R, load the `NHANES` package, and then run the command `data(NHANES)` which will load the NHANES data. Type `?NHANES` and read about the dataset.
```{r}
# install.packages("NHANES")
library(NHANES)
data(NHANES)
?NHANES
```

2. Make an object `nhanes` that is a subset version `NHANESraw` that does not include any missing data for `Diabetes`, `BPSysAve`, `BPDiaAve`, or `Age`.
```{r}
library(tidyr)
library(dplyr)
nhanes <- NHANESraw %>% filter(!is.na(Diabetes) & !is.na(BPSysAve) & !is.na(BPDiaAve) & !is.na(Age))
```

3. (1 point) Further subset the data such the observations with `BPDiaAve` equal to zero are removed.
```{r}
nhanesnozeros <- nhanes %>% filter(BPDiaAve != 0)
```

4. (1 point)
Make an object `nhanes09` that is a subset of `nhanes` to only the 2009_10 data. This will be your training dataset. Also make an object `nhanes11` that is a subset of `nhanes` to only the 2011_12 data. This will be your test dataset.
```{r}
nhanes09 <- nhanesnozeros %>% filter(SurveyYr == "2009_10")
nhanes11 <- nhanesnozeros %>% filter(SurveyYr == "2011_12")
#rm(nhanesnozeros)  # to save memory
#rm(nhanes)  # to save memory
```

## Logistic regression

5. (2 point) Fit a logistic regression model (call it `glm1`) using the `nhanes09` dataset. Use `Diabetes` as the outcome and averaged systolic blood pressure (`BPSysAve`) as a single predictor. Use the summary command to examine the fitted model. Generate the 95% confidence intervals for the `BPSysAve` coefficient.
```{r}
glm1 <- glm(Diabetes ~ BPSysAve, family="binomial", data=nhanes09)
summary(glm1)
confint(glm1)
```

6. (1 point) Generate the estimate and 95% confidence interval for the odds-ratio associated with BPSysAve. Summarize the result.
```{r}
exp(coefficients(glm1)["BPSysAve"])
exp(confint(glm1)["BPSysAve",])

# For each unit increase in Systolic Blood Pressure we see about a 2.93% increase in odds of diabetes, 95% CI 2.58% to 3.29%.
```

7. (1 point) Predict the probabilities of diabetes associated with each of the training observations of `BPSysAve`. Make a vector of predictions for diabetes based on whether the predictions are above or below 0.5.
```{r}
pred09 <- predict(glm1, type="response")  ## Generate predicted probabilities
pred09binary <- ifelse(pred09 > 0.5, "Yes", "No") ## Threshold at 0.5
```

8. (1 point) Generate a confusion matrix that shows the number of false positives, false negatives, true positives, and true negatives in the training data. The rows should correspond to the true diabetes status and the columns should correspond to the predicted values.
```{r}
table(true_disease=nhanes09$Diabetes, predictions=pred09binary)
```

9. (1 point) Find the proportion of correctly classified observations in the training data.
```{r}
mean(nhanes09$Diabetes == pred09binary)
```

10. (2 points) Now repeat questions 7 to 9 but for predicting the test dataset.
```{r}
pred11 <- predict(glm1,type="response",newdata=nhanes11)
pred11binary <- ifelse(pred11 > 0.5,"Yes","No")
xt11 <- table(true_disease=nhanes11$Diabetes, predictions=pred11binary)
print(xt11)
print(paste("Test accuracy", mean(nhanes11$Diabetes == pred11binary)))
```

11. (1 point) Comment on the difference in results between the training and test prediction tables and classification accuracies.
```{r}
# The results are actually very similar for the training and test datasets. This is not totally surprising given that nearly all predicted probabilities are below 0.5.
```

12. (1 point) Manually calculate the sensitivity and specificity estimates for the test dataset based on the 0.5 threshold.
```{r}
sens <- xt11[2,2]/sum(xt11[2, ])
spec <- xt11[1,1]/sum(xt11[1, ])
print(paste("Sensitivity", sens))
print(paste("Specificity", spec))
```

13. (2 points) Generate an ROC curve using the test data. What is the AUC and its 95% confidence interval?
```{r}
library(ROCit)
ROCit_obj <- rocit(score=pred11, class=nhanes11$Diabetes)
plot(ROCit_obj)
ciAUC(ROCit_obj)
```

14. What value can you use to threshold the predicted probability to achieve a sensitivity of at least 0.6 and a specificity of at least 0.3?
```{r}
# The following cutoffs can all be used to achieve the desired sensitivity and specificity levels
thresholds <- ROCit_obj$Cutoff[(ROCit_obj$TPR > 0.6) & (ROCit_obj$FPR < 0.3)]
print(paste("Possible threshold:", thresholds[1]))
```

15. (2 points) Comment on the results of the analyses for the different thresholds in terms of the tables, classification accuracies, and sensitivity and specificity. Under what circumstances might you prefer each of the thresholds?

```{r}
# The accuracy is highest for the 0.5 threshold, is a little worse for the 0.2 threshold but quite a bit worse for the 0.1 threshold. However, for the 0.5 and 0.2 thresholds the specificity is really high but the sensitivity is poor (slightly better sensitivity for 0.2). Despite having considerably lower classification accuracy, there is more of a balance between sensitivity and specificity for the 0.1 threshold. This set of results reflects the fact that a large majority of the subjects do not have diabetes so in the absence of strong evidence for diabetes you are more likely to get it right by predicting no diabetes. However, if you want to screen for diabetes you need to get the sensitivity up and so may want to use the 0.1 threshold (albeit at the risk of having many false positives.)
```

16. (2 points) Fit a multiple predictor logistic regression (call it `glm2`) with `Diabetes` as outcome and predictors: `BPSysAve`, `BPDiaAve`, and `Age`. Use the `summary` command to examine the fitted model and determine the estimated coefficients, odds-ratios, and 95% confidence intervals thereof.
```{r}
glm2 <- glm(Diabetes ~ BPSysAve + BPDiaAve + Age, family="binomial", data=nhanes09)
summary(glm2)
# The 95% confidence interval for the model coefficients
confint(glm2)
# The odds-ratios associated with each predictor
exp(coefficients(glm2))
# The 95% confidence interval for the odds-ratios associated with each predictor
exp(confint(glm2))
```

17. (2 points) Generate an ROC curve for the `glm2` model using the test data. What is the AUC and its 95% confidence interval?
```{r}
pred11_glm2 <- predict(glm2, newdata = nhanes11, type = "response")
ROCit_obj2 <- rocit(score=pred11_glm2, class=nhanes11$Diabetes)
plot(ROCit_obj2)
ciAUC(ROCit_obj2)
```

18. (1 point) What is the maximum sensitivity level you can achieve if we require the specificity to be at least 0.3?
```{r}
max(ROCit_obj2$TPR[ROCit_obj2$FPR < 0.7])
```

19. (1 point) Would you prefer the single predictor or multiple predictor model if your objective was to maximize classification accuracy, and which threshold level would you choose? Comment on the reason for your choices.

```{r}
# Probably would threshold based on 0.5 since that is the best threshold for classification accuracy for both models. You might prefer the multi predictor model since it has very slightly higher accuracy on the test set, but it is so marginal that you may opt for the simpler single predictor model, even with the tiny higher performance on the test set.
```

## Linear discriminant analysis

20. (2 points) Fit a linear discriminant analysis (`lda1`) with `Diabetes` as outcome and predictors of `BPSysAve`, `BPDiaAve`, and `Age` in the training dataset. Examine the fit by typing `lda1`.
```{r}
lda1 <- lda(Diabetes ~ BPSysAve + BPDiaAve + Age, data=nhanes09)
lda1
```

21. (2 points) Generate the confusion matrix for `lda1` using the test set. Compute the classification accuracy, sensitivity, and specificity.
```{r}
predlda1 <- predict(lda1, newdata=nhanes11)
lda1tab <- table(true_disease=nhanes11$Diabetes, predictions=predlda1$class)
print(lda1tab)
print(paste("Accuracy", mean(predlda1$class==nhanes11$Diabetes)))
senslda1 <- lda1tab[2,2]/sum(lda1tab[2, ])
speclda1 <- lda1tab[1,1]/sum(lda1tab[1, ])
print(paste("Sensitivity", senslda1))
print(paste("Specificity", speclda1))
```

22. How do these measures compare with that of the logistic regression model with these predictors and 0.5 threshold?
```{r}
# The accuracy is pretty similar 0.8869 (lda) vs. 0.8895 (logistic) between the two approaches. The lda has higher sensitivity relative to logistic at the expense of some specificity.
```


23. (3 points) Redo question 21 but with prior probabilities set to 0.5 for diabetes.
```{r}
predlda2 <- predict(lda1, prior=c(0.5,0.5), newdata=nhanes11)
lda2tab <- table(true_disease=nhanes11$Diabetes, predictions=predlda2$class)
print(lda2tab)
print(paste("Accuracy", mean(predlda2$class==nhanes11$Diabetes)))
senslda2 <- lda2tab[2,2]/sum(lda2tab[2, ])
speclda2 <- lda2tab[1,1]/sum(lda2tab[1, ])
print(paste("Sensitivity", senslda2))
print(paste("Specificity", speclda2))
```

24. (2 points) Comment on how LDA's performance changed when we changed the prior probabilities.
```{r}
# The second model assumedtest patients had an equal probability of being diabetic and non-diabetic a priori. This led to an increase in sensitivity but decreased the specificity. The overall accuracy decreased because the prior distribution did not reflect the proportions seen in the test data.
```


## Penalized regression

26. Read in the dementia data "dementia2.csv" into a data frame called `dementia_dat`. This dataset contains measurements obtained from MRI brain scans and whether or not the patient has dementia. We'll try to build a prediction model for diagnosing dementia based on these derived measurements. How many observations are in this dataset? How many predictors are in this dataset?
```{r}
dementia_dat <- read.csv("dementia2.csv")
print(paste("Number of observations", nrow(dementia_dat)))
print(paste("Number of predictors", ncol(dementia_dat) - 1))
head(dementia_dat)
```

27. Load the `glmnet` and `caret` packages.
```{r}
library(glmnet)
library(caret)
```

28. (1 point) Set the random seed to 4 and then split the data into 2 sets (400 train, and 260 test) 
```{r}
set.seed(4)
groups <- c(rep("train",400),rep("test",260))
groups <- sample(groups,length(groups))
datTrain <- dementia_dat[groups=="train", ]
datTest <- dementia_dat[groups=="test", ]
```

29. (4 points) Perform cross-validated lasso in the training data to select the optimal penalty parameter lambda. Use 5 folds and search over the range $\lambda = 10^{3}$ to $\lambda = 10^{-3}$. Set `Dementia` as outcome with all other variables as predictors except `MacCohort_kr`. Use the `caret` package to do CV.
```{r}
train_control <- trainControl(method="cv", number=5)
grid <- 10^seq(3,-3,length=100)
caret_grid <- data.frame("lambda" = grid, "alpha"= 1)
cv_model <- train(Dementia ~ . - MacCohort_kr, data=datTrain, trControl=train_control, method="glmnet", tuneGrid=caret_grid)
```

30. (1 point) What is the optimal value of lambda?
```{r}
cv_model$bestTune$lambda
```


31. (1 point) Generate the confusion matrix for this final model.
```{r}
lassoPred <- predict(cv_model, newdata = datTest)
table(true_disease=datTest$Dementia, predictions=lassoPred)
```

32. (1 point) How many non-zero coefficients are in the final model?
```{r}
lasso_coefs <- coefficients(cv_model$finalModel, s = cv_model$bestTune$lambda)
sum(abs(lasso_coefs) > 0)
```