HSML 6295 s3 Software Practice ANSWERS.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### HSML 6295 Session 3 (Classification) - ANSWERS\n",
                "\n",
                "#### I. The \"Framingham Heart Study\" Data Set\n",
                "\n",
                "\"Cardiovascular disease (CVD) is the leading cause of death and serious illness in the United States. In 1948, the Framingham Heart Study - under the direction of the National Heart Institute (now known as the National Heart, Lung, and Blood Institute or NHLBI) - embarked on an ambitious project in health research. ...\n",
                "The researchers recruited 5,209 men and women between the ages of 30 and 62 from the town of Framingham, Massachusetts, and began the first round of extensive physical examinations and lifestyle interviews that they would later analyze for common patterns related to CVD development. Since 1948, the subjects have continued to return to the study every two years for a detailed medical history, physical examination, and laboratory tests.\" \n",
                "https://framinghamheartstudy.org/fhs-about/history/\n",
                "\n",
                "Our objective is to predict the 10-year risk of Coronary Heart Disease (CHD). The response variable `observed.CHD` is 1 if the individual developed CHD over the 10-year observation period and 0 otherwise.\n",
                "\n",
                "| Predictor | Description\n",
                "| ---       | ---\n",
                "| Male | sex of patient\n",
                "| Age | age in years at first examination\n",
                "| Education | some HS [1], high school [2], some college [3], college [4]\n",
                "| Smoker | yes [1], no [0]\n",
                "| Cigs.per.Day | cigarettes smoked per day\n",
                "| BP.Medication | on blood pressure medication at time of first examination\n",
                "| Stroke | has had a stroke\n",
                "| Hypertension | currently hypertensive\n",
                "| Diabetes | currently has diabetes \n",
                "| Cholesterol | total cholesterol (mg/dL)\n",
                "| Systolic.BP | systolic blood pressure\n",
                "| Diastolic.BP | diastolic blood pressure\n",
                "| BMI | Body Mass Index, weight (kg)/height (m)^2\n",
                "| Heart.Rate | heart rate (beats/minute)\n",
                "| Glucose | blood glucose level (mg/dL)\n",
                "\n",
                "Read in the data set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "framingham = read.csv(\"HSML 6295 ds Framingham.csv\")\n",
                "dim(framingham)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "The `na.omit()` function removes all the observations with missing values for any variable.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "framingham = na.omit(framingham)\n",
                "dim(framingham)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "[Load the R package \"stargazer\" for displaying the summary statistics succinctly:]\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(stargazer)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "stargazer(framingham, \n",
                "          summary.stat = c(\"n\", \"mean\", \"sd\", \"min\", \"p25\", \"median\", \"p75\", \"max\"),\n",
                "          type = \"text\", title=\"Descriptive statistics\", digits=1)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Tabulate relative frequencies of the values of `observed.CHD` in the full data set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "round(prop.table(table(framingham$observed.CHD)),2)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Declare `observed.CHD` a factor variable and assign labels \"Yes\" and \"No\" to 1 and 0, respectively:\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "framingham$observed.CHD = factor(framingham$observed.CHD, \n",
                "                             levels = c(0,1), labels = c(\"No\", \"Yes\"))\n",
                "round(prop.table(table(framingham$observed.CHD)),2)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Randomly split the `framingham` data set into a *training set* comprising 75% of observations and a *test set* comprising the remainder.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(caTools)\n",
                "set.seed(1)\n",
                "split = sample.split(framingham$observed.CHD, SplitRatio = 0.75)\n",
                "train = subset(framingham, split==TRUE)\n",
                "test = subset(framingham, split==FALSE)\n",
                "dim(framingham)\n",
                "dim(train)\n",
                "dim(test)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### II. Baseline Prediction \n",
                "\n",
                "Without information about any predictors, the best approach is to classify all observations to the most common class in the training data.\n",
                "\n",
                "Tabulate relative frequencies of the values of `observed.CHD` in the training set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "round(prop.table(table(train$observed.CHD)),2)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 1.** Declare `Education` a factor variable and assign labels \"some HS\", \"high school\", \"some college\", \"college\" to values 1, 2, 3, and 4, respectively, in the training set. Then tabulate relative frequencies of the values of `Education` in the training set. What is the baseline prediction for `Education`?\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "train$Education = factor(train$Education,     # declare Education as factor \n",
                "                         levels = c(1, 2, 3, 4), \n",
                "                         labels = c(\"some HS\", \"high school\", \"some college\", \"college\"))\n",
                "round(prop.table(table(train$Education)),2)   # tabulate relative frequencies\n",
                "train$Education = as.numeric(train$Education) # convert Education from factor back to numeric format\n",
                "round(prop.table(table(train$Education)),2)   # tabulate relative frequencies\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "    \n",
                "**The baseline prediction for `Education` is \"some HS\" because it is the most common class in the training set.**\n",
                "\n",
                "Create a vector of baseline (\"naïve\") predictions: \n",
                "\n",
                "Every individual in the test set is predicted to have \"No\" CHD.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "predicted.CHD = rep(\"No\",nrow(test))\n",
                "round(prop.table(table(predicted.CHD)),2)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Confusion matrix\n",
                "\n",
                "Refer to Table 4.6 in ISLR.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "|                 |         |                |       |       |\n",
                "| ---             | ------- | ---:           | ---:  | ---:  |\n",
                "|                 |         | **Predicted**  | **CHD**  \n",
                "|                 |         | No             | Yes   | Total |\n",
                "|**Observed CHD** | No      | 775 TN         | 0 FP  | 775 N |\n",
                "|                 | Yes     | 139 FN         | 0 TP  | 139 P |\n",
                "|                 | Total   | 914 N*         | 0 P*  | 914 A |\n",
                "\n",
                "\n",
                "We compute 3 statistics that measure predictive performance\n",
                "\n",
                "Statistic  \n",
                "Definition  \n",
                "Value  \n",
                "\n",
                "**Accuracy**  (1 - Error Rate)  \n",
                "(TN+TP)/A  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+0)/cm[3,2],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**True Positive Rate**  (Sensitivity)  \n",
                "TP/P  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*0/cm[2,2],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**False Positive Rate**  (1 - Specificity)  \n",
                "FP/N  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*0/cm[1,2],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Note**  \n",
                "Here's a mnemonic that I use to remember the definitions of sensitivity and specificity:  \n",
                "The word \"sensitivity\" includes the letter \"n\" but is defined as the true **p**ositive rate.  \n",
                "Conversely, the word \"specificity\" includes the letter \"p\" but is defined as the true **n**egative rate.\n",
                "\n",
                "\n",
                "#### III. Logistic Regression\n",
                "\n",
                "We obtain the predictions in 3 steps:\n",
                "\n",
                "1. Fit the logistic regression model using the observations in the training set\n",
                "2. Use the fitted model to estimate the probability of CHD for each observation in the test set\n",
                "3. Use the estimated probabilities and a threshold value to classify the test observations\n",
                "\n",
                "**Step 1.** Estimate a logistic regression model of `observed.CHD` as a function of all predictors using the observations in the *training* set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fit = glm(observed.CHD ~ ., data = train, family=binomial)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Print the results\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "summary(fit)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**Step 2.** Estimate the probability of CHD $\\hat{p}$ for each observation in the *test* set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "probability.CHD = predict(fit, type=\"response\", newdata=test)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**Step 3.** Classify each observation in the *test* set for which the probability is greater than 0.5, i.e. $\\hat{p} > 0.5$, as `predicted.CHD = \"Yes\"`; else classify the observation as `predicted.CHD = \"No\"`. \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.5] = \"Yes\"\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "The predictions from Step 3 yield the confusion matrix and predictive performance statistics\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE) # add row & column totals\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Definition  \n",
                "Value  \n",
                "\n",
                "**Accuracy**  \n",
                "(TN+TP)/A  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**True Positive Rate**  \n",
                "TP/P  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**False Positive Rate**  \n",
                "FP/N  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 2.** Compute the confusion matrix and the 3 resulting predictive performance statistics (accuracy, true positive rate, false positive rate) when the probability threshold is 0.48.\n",
                "\n",
                "**The confusion matrix with threshold of 0.48 and the resulting predictive performance statistics are**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.48] = \"Yes\"\n",
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**  \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 3.** Can you find a probability threshold for which the true positive rate exceeds 50%. What is the classifier's false positive rate in this case?\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.20] = \"Yes\"\n",
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**A threshold of 0.20 achieves a true positive rate of**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "**but also raises the false positive rate to**\n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Varying the probability threshold yields the ROC curve and AUC (area under the ROC curve)\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc_full = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE, \n",
                "               xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "               auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Note**  \n",
                "Here's a mnemonic that I use to remember which rates are shown on the horizontal and vertical axes, respectively:  \n",
                "The first letter (F) in \"False Positive Rate\" includes includes two horizontal lines. We show the FPR on the horizontal axis.  \n",
                "The first letter (T) in \"True Positive Rate\" is composed of two lines, of which the vertical line is longer than the horizontal line. We show the TPR on the vertical axis.\n",
                "\n",
                "\n",
                "**Knowledge Check 4.** Estimate a logistic regression model on the training data but only include the following 5 predictors: Male, Age, Cigs.per.Day, Systolic.BP, Glucose. These are the predictors whose p-values were smaller than 0.05 (had at least one asterisk) in the results table above. Compute the confusion matrix and the 3 resulting predictive performance statistics (accuracy, true positive rate, false positive rate) when the probability threshold is 0.48, 0.47, 0.46, 0.45, and 0.44. If missing a case of CHD was much more costly than mistakenly diagnosing a patient with CHD, which threshold value would you choose?\n",
                "\n",
                "**Run the following chunk and only replace the threshold value in the following line of code:**\n",
                "`predicted.CHD[probability.CHD > 0.44] = \"Yes\"`\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fit = glm(observed.CHD ~ Male + Age + Cigs.per.Day + Systolic.BP + Glucose, \n",
                "          data = train, family=binomial)\n",
                "probability.CHD = predict(fit, type=\"response\", newdata=test)\n",
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.44] = \"Yes\"\n",
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n",
                "accuracy = round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "accuracy\n",
                "tpr = round(100*cm[2,2]/cm[2,3],1)\n",
                "tpr\n",
                "fpr = round(100*cm[1,2]/cm[1,3],1)\n",
                "fpr\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**For all 5 threshold values, the accuracy is the same. As we lower the threshold, we predict more CHD cases. Thus, we raise both the true positive rate and the false positive rate. If raising the true positive rate is more valuable than lowering the false positive rate, we would prefer the lowest threshold that still maintains the same accuracy level, i.e. 0.44.**\n",
                "\n",
                "**Knowledge Check 5.** Construct the ROC curve and compute the AUC. How does the AUC of the restricted model (the one using only 5 predictors) compare to the AUC of the full model? Now change the value of the seed in line 100 from 1 to 1000 and compare again the AUC values for the full and restricted models.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc_restricted = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE, \n",
                "                     xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "                     auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**When the seed is 1, the AUC of the full model is**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(roc_full$auc, 3)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "**while the AUC of the restricted model is** \n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(roc_restricted$auc, 3)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**When the seed is 1000, the AUC of the full model is 0.734, while the AUC of the restricted model is 0.739.**\n",
                "\n",
                "#### IV. Linear Discriminant Analysis (LDA)\n",
                "\n",
                "As with logistic regression, with LDA we obtain the predictions in 3 steps:\n",
                "\n",
                "1. Fit the model on the training set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(MASS)\n",
                "fit = lda(observed.CHD ~ ., data = train)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "2. Compute the estimated probability of CHD for each observation in the test set\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "probability.CHD = predict(fit, newdata = test)$posterior[,2]\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "3. Use the estimated probabilities and a threshold value to classify the test observations\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.5] = \"Yes\"\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "The predictions from Step 3 yield the confusion matrix and predictive performance statistics\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Varying the probability threshold yields the ROC curve and AUC (area under the ROC curve)\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE, \n",
                "          xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "          auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 6.** Rerun the linear discriminant analysis on the training data but only include the following 5 predictors: Male, Age, Cigs.per.Day, Systolic.BP, Glucose. Compute the confusion matrix and the 3 resulting predictive performance statistics (accuracy, true positive rate, false positive rate) when the probability threshold is 0.5. Construct the ROC curve and compute the AUC.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fit = lda(observed.CHD ~ Male + Age + Cigs.per.Day + Systolic.BP + Glucose, \n",
                "          data = train)\n",
                "probability.CHD = predict(fit, newdata = test)$posterior[,2]\n",
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.5] = \"Yes\"\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE, \n",
                "          xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "          auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### V. Quadratic Discriminant Analysis (QDA)\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(MASS)\n",
                "fit = qda(observed.CHD ~ ., data = train)\n",
                "probability.CHD = predict(fit, newdata = test)$posterior[,2]\n",
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.5] = \"Yes\"\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE, \n",
                "          xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "          auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 7.** Rerun the quadratic discriminant analysis on the training data but only include the following 5 predictors: Male, Age, Cigs.per.Day, Systolic.BP, Glucose. Compute the confusion matrix and the 3 resulting predictive performance statistics (accuracy, true positive rate, false positive rate) when the probability threshold is 0.5. Construct the ROC curve and compute the AUC.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fit = qda(observed.CHD ~ Male + Age + Cigs.per.Day + Systolic.BP + Glucose, \n",
                "          data = train)\n",
                "probability.CHD = predict(fit, newdata = test)$posterior[,2]\n",
                "predicted.CHD = rep(\"No\", nrow(test))\n",
                "predicted.CHD[probability.CHD > 0.5] = \"Yes\"\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "cm <- addmargins(cm, FUN = list(Total = sum), quiet = TRUE)\n",
                "cm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Statistic  \n",
                "Value\n",
                "\n",
                "\n",
                "**Accuracy**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*(cm[1,1]+cm[2,2])/cm[3,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**True Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[2,2]/cm[2,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "**False Positive Rate**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                " round(100*cm[1,2]/cm[1,3],1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(pROC)\n",
                "roc = roc(test$observed.CHD, probability.CHD, plot=TRUE, legacy.axes = TRUE,\n",
                "          xlab = \"False Positive Rate\", ylab = \"True Positive Rate\",\n",
                "          auc.polygon=TRUE, max.auc.polygon=TRUE, print.auc=TRUE)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### VI. K-Nearest Neighbors\n",
                "\n",
                "Convert the factor variable `observed.CHD` back to numeric format:\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "train$observed.CHD = as.numeric(train$observed.CHD)-1\n",
                "test$observed.CHD = as.numeric(test$observed.CHD)-1\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Define a function that computes the accuracy statistic\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy = function(Actual, Predicted) {\n",
                "    round(100*mean(Actual == Predicted),2)\n",
                "}\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "Run K-NN and print the resulting confusion matrix and accuracy value\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(class)\n",
                "predicted.CHD = knn(train[,-1], test[,-1], train[,ncol(train)], k = 1)\n",
                "cm = table(test$observed.CHD, predicted.CHD)\n",
                "(cm = addmargins(cm, FUN = list(Total = sum), quiet = TRUE))\n",
                "(accuracy_1 = accuracy(Actual = test$observed.CHD, Predicted = predicted.CHD))\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Knowledge Check 8.** Compute the accuracy of the K-nearest neighbor classifier for K = 2, 4, 8, 16, 32, 64, and 128 and record it in the table below. Which value of K minimizes the test error rate?\n",
                "\n",
                "K (number of nearest neighbors used in prediction)\n",
                "Accuracy\n",
                "\n",
                "1\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_1\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Replace the value of `k` in the line of code**\n",
                "\n",
                "`predicted.CHD = knn(train[,-1], test[,-1], train[,ncol(train)], k = 1)`\n",
                "\n",
                "**above with 2, 4, 8, 16, 32, 64, and 128:**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "library(class)\n",
                "k_values = list(1, 2, 4, 8, 16, 32, 64, 128)\n",
                "for(i in k_values){\n",
                "  predicted.CHD = knn(train[,-1], test[,-1], train[,ncol(train)], k = i)\n",
                "  (accuracy_k = accuracy(Actual = test$observed.CHD, Predicted = predicted.CHD))\n",
                "  result = paste(\"accuracy_\", i, sep = \"\")\n",
                "  assign(result, accuracy_k)\n",
                "  }\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "K (number of nearest neighbors used in prediction)\n",
                "Accuracy\n",
                "\n",
                "1\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_1\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "2\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_2\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "4\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_4\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "8\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_8\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "16\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_16\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "32\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_32\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "64\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_64\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "\n",
                "128\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "accuracy_128\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": "\n"
        }
    ],
    "metadata": {
        "anaconda-cloud": "",
        "kernelspec": {
            "display_name": "R",
            "langauge": "R",
            "name": "ir"
        },
        "language_info": {
            "codemirror_mode": "r",
            "file_extension": ".r",
            "mimetype": "text/x-r-source",
            "name": "R",
            "pygments_lexer": "r",
            "version": "3.4.1"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 1
}