---
title: "Take-Home Final Exam"
author: "Statistical Computing, 36-350"
date: "Monday December 3, 2018"
---

```{r, include=FALSE}
knitr::opts_chunk$set(cache=TRUE, autodep=TRUE, cache.comments=TRUE)
```

Name:  
Andrew ID: 

Instructions
===

- You must work alone.

- You will be submitting a knitted HTML file to Canvas, as usual, by **Sunday December 9** at 11:59pm.

- For full credit on each question, make sure you follow exactly what is asked, and answer each prompt. The total is 66 points (and there are 10 bonus points).

- You may only ask clarification questions on Piazza; you may not ask questions for help. This is a final exam, not a lab assignment.

- This should go without saying (but we have had several problems in past years): do not cheat. It ends poorly for everyone.

Simulation and plotting (18 points)
===

- **1a (4 points).** Simulate 1000 pairs $(X_1,X_2)$ according to
$$
X_1 \sim N(0,1.5) \quad \text{and} \quad X_2 \sim N(0,1), \quad \text{independently}.
$$
(Hhere $N(\mu,\sigma^2)$ denotes the normal distribution with mean $\mu$ and variance $\sigma^2$.) Call this data "group A", and save the data points in matrix called `a.mat`, of dimension 1000 x 2. Simulate 1000 pairs $(X_1,X_2)$ according to
$$
X_1 \sim N(2.5,1) \quad \text{and} \quad X_2 \sim N(3,1.5), \quad \text{independently}.
$$
Call this data "group B", and save the data points in a matrix `b.mat`, of dimension 1000 x 2. Compute the column means and column standard deviations of `a.mat` and `b.mat`, and show that they're close to what you would expect.

- **1b (4 points).** On a single plot, plot $X_2$ versus $X_1$, using all the points from groups A and B. Color the points from group A in black, and those from group B in red. Set the x label to "X1", the y label to "X2", and leave the title blank. Make sure the axes limits are just big enough to show all the points on the plot, and set these programmatically (do not set these using hard constants). Add a legend to explain the color coding. Do you see significant overlap in the points from the two groups?

- **1c (3 points).** Define a data frame `dat` of dimension 2000 x 3, as follows. First, stack the rows of `a.mat` and `b.mat` together, to get a data frame of dimension 2000 x 2. Then, append a column of 0s and 1s, that indicates whether the points (rows) come from group A or B (with 0 for A and 1 for B). Finally, name the columns of `dat` to be `X1`, `X2`, `Y` respectively. Show the first 3 and last 3 rows of `dat`.

- **1d (2 points).** Perform a logistic regression of `Y` (response) on `X1` and `X2` (variables), using the data `dat`. Use the function `glm()` with `family="binomial"`, and save the output object as `log.mod`. Allow for an intercept in the model (the default). Report the intercept, and the coefficients for `X1` and `X2`.

- **1e (5 points).** Let $\beta_1,\beta_2$ be the coefficients of $X_1,X_2$ that you computed in the previous logisitic regression, and let $\beta_0$ be the intercept. Consider the line defined by the equation 
$$
\beta_0 + \beta_1 X_1 + \beta_2 X_2 = 0.
$$
Via simple algebra, rewrite this equation as 
$$
X_2 = a + b X_1,
$$
for a particular $a,b$ that you will compute. Then reproduce your plot from Q1b, and plot the line $X_2 = a + b X_1$ on top, as a thick blue dashed line. (In case you're curious, this is called the *decision boundary* from your logistic regression.) Modify your legend so that your line explained. Does your line look like a reasonable separator of groups A and B?

Iteration (14 points)
===

- **2a (2 points).** Refer back to the logistic regression object you computed in Q1d. Using `summary()`, form a vector of the standard errors of the intercept and coefficients. Display this vector (it should be of length 3).

- **2b (10 points).** Now we're going perform an iterative strategy to estimate these standard errors ourselves. It works as follows:
    i. Sample 2000 $(X_1,X_2,Y)$ triplets with replacement from your original data set of 2000 $(X_1,X_2,Y)$ triplets (the data set with groups A and B together).
    ii. Run a logistic regression of $Y$ on $X_1,X_2$ (including intercept), using the data constructed in from Step i. Store the coefficient vector.
    iii. Repeat Steps i and ii 100 times. You will have 100 coefficient vectors. Report the standard deviation of the intercept coefficients, and the $X_1$ coefficients, and the $X_2$ coefficients.

    Perform this iterative strategy (in case you're curious, this is known as *bootstrapping*) on your data set from Q1c. Hint: a `for()` loop is probably simplest. What are the standard deviations you get for the intercept and $X_1$ and $X_2$ coefficients? Are they close to the standard errors you found in Q2a? (They should be.)

- **2c (2 points).** Rerun your code from Q2b, but now sampling *without* replacement (before you were sampling with replacement). What are the standard errors you estimate in the end? And does this make sense to you?

Functions and debugging (16 points)
===

- **3a (3 points).** Below is a function `predict.fun()` that your instructor wrote to get predictions out of a logistic regression model that was fit on two variables. The user passes in `log`, a fitted logistic regression object (from `glm()`), and `x1` and `x2`, values of the two variables at which he/she wants to make predictions. The function is supposed to return the predicted class (0 or 1), at each value of `x1` and `x2` that is passed. Note that `x1` and `x2` can be vectors; so, e.g., if these are of length 5, then the function should also return a vector of length 5, giving the class predictions for each of the 5 pairs of `x1` and `x2` points.

    Unfortunately, your instructor was sloppy and introduced bugs in his implementation. Fix them, then check that the function returns the outputs as specified in the comments for the examples below. Note that `log.mod` is your fitted logistic regression object from Q1d. When done, remove `eval=FALSE` in the option for the code chunk. Hint: the argument `type="response"` in the call to `predict` is *not* a bug. This is making sure that the `predict()` function returns probabilities. These are are from the logistic model; if the probability is above 1/2, we predict class 1, and if it is below 1/2, we predict class 0.

```{r, eval=FALSE}
predict.fun = function(mod, x1, x2) {
  prob = predict(mod, newdata=data.frame(x1=x1, x2=x2), type="response")
  if (prob <= 0.5) return(0)
  else return(1)
}

predict.fun(log.mod, -10, -10) # Should return 0
predict.fun(log.mod, 10, 10) # Should return 1
predict.fun(log.mod, c(-10,10), c(-10,10)) # Should return c(0,1)
```

- **3b (5 points).** Below is an alternative function `predict.fun2()` that your instructor wrote that serves the same purpose (and has the same argument structure) as `predict.fun()`, but uses a different strategy. The previous one calls `predict()` on a fitted logistic regression object; this one directly extracts the coefficients from the logistic regression object, call them $\beta_0,\beta_1,\beta_2$ in math notation, and then computes the linear combination
$$
\beta_0 + \beta_1 X_1 + \beta_2 X_2
$$
for each value of the variables $X_1$ and $X_2$ in consideration. It then compares this linear combination to 0: if below 0, it predicts class 0, and if above 0, it predicts class 1. There are several bugs in `predict.fun2()`. Fix them, check that you get the outputs as specified in the comments for the examples below, then remove `eval=FALSE`.

```{r, eval=FALSE}
predict.fun2 = function(mod, x1, x2) {
  link = sum(c(x1,x2) * coef(mod))
  if (link <= 0) return(0)
  else return(1)
}

predict.fun2(log.mod, -10, -10) # Should return 0
predict.fun2(log.mod, 10, 10) # Should return 1
predict.fun2(log.mod, c(-10,10), c(-10,10)) # Should return c(0,1)
```

- **3c (8 points).** Below is a function `plot.predictions()` that your instructor wrote that is supposed to take a fitted logistic regression model `mod`, and make predictions from this model over a grid of `x1` and `x2` values, where the grid is 100 x 100, and has endpoints defined by `x1.lim` and `x2.lim`. The function is then supposed to produce a plot that colors the points on this grid by their predicted classes, with `col0` for class 0 and `col1` for class 1 (these are both strings to be passed in by the user, but with default values of "black" and "red", respectively).

    Your instructor didn't introduce any bugs (phew!) but stopped short and didn't finish the implementation. Complete the missing code, as indicated by the comments. Then run it with appropriate arguments to produce a plot based on your logistic regression model `log.mod` from Q1d. Set the grid endpoints `x1.lim` and `x2.lim` so that your plot is just big enough to capture all the data from Q1a. Then, on top of the plot you get from `plot.predictions()`, plot the data points themselves, and the decision boundary, just as you did in Q1e (you may copy over your code from Q1e here; just make sure everything is on one plot in the end). 
    
```{r}
plot.predictions = function(mod, x1.lim, x2.lim, col0="black", col1="red") {
  # Define x1 and x2 variables over the grid. Fill in the ?? below, then 
  # uncomment the lines
  #x1.grid = seq(??, ??, len=100)
  #x2.grid = seq(??, ??, len=100)
  x1 = rep(x1.grid, times=100)
  x2 = rep(x2.grid, each=100)
  
  # Call predict.fun() to get the predicted classes. (You can change this
  # to predict.fun2() if you want; just use whichever one you got working)
  pred = predict.fun(mod, x1, x2)
  
  # Define the colors based on the predicted classes. Remember that col0 
  # is for class 0, and col1 is for class 1. Fill in the ??, uncomment
  #cols = ifelse(??, ??, ??)
  
  # Plot the grid points with the right colors. Fill in the ??, uncomment
  #plot(??, ??, col=cols, xlab="X1", ylab="X2", pch=20, cex=0.1)
}
```
    
Data cleaning and tidyverse (18 points)
===

- **4a (2 points).** A data set on horse colic is available at https://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.data. Read this into your R session using `read.table()`; make sure to set `stringsAsFactors=FALSE`; save the resulting data frame as `horse.df`; check that it has dimension 300 x 28; and set its column names to the string vector `nms` below. Print out the first 5 rows and (all columns) of `horse.df` to the console. To read more about the data set, visit https://archive.ics.uci.edu/ml/datasets/Horse+Colic.

```{r}
nms = c("Surgery", "Age", "Hosp_ID", "Temp_Rect", "Pulse", "Resp", "Temp_Extr",
        "Pulse_Peri", "Mucous", "Cap_Time", "Pain", "Peris", "Ab_Dist", "Tube",
        "Reflux", "Reflux_PH", "Feces", "Abdomen", "Cell_Vol", "Total_Protein",
        "Ab_Appear", "Ab_Total_Prot", "Outcome", "Lesion", "Lesion_Site", 
        "Lesion_Type", "Lesion_Subt", "Lesion_Code")
```
 
- **4b (3 points).** How many columns in `horse.df` are of character type, and of integer type? How many columns have missing values (represented by "?" in this data set), and what are they (what are their columns names)? The UCI website reports that about 30% of this data set is missing. Is this accurate? What do you calculate the missingness percentage to be?

- **4c (3 points).** Convert all of the "?" values in `horse.df` to NA values. After this, convert all of the character columns to numeric, but make sure that the `horse.df` object is still a data frame (just explicitly cast it to one, if you need to; otherwise you won't be able to use tidyverse + pipes in what follows). When done, print the first 5 rows (and all columns).

- In the following questions, use the tidyverse (specifically, `dplyr` verbs) along with pipes to construct your answer. Each part should require just a single flow of pipes.

- **4d (3 points).** Arrange the horses (the rows in `horse.df`) by `Surgery` (1 means yes surgery, 2 means no surgery), and then by descending values of `Pulse`. Report the `Surgery`, `Pulse`, `Temp_Rect` (rectal temperature), and `Ab_Dist` (abdominal distention: 1 being best, 4 being worst) variables, for all horses with `Pulse` at least 130.

- **4e (3 points).** For each `Pain` level (1 being best, 5 being worst), compute the average `Pulse` (and make sure to exclude NA values in this calculation), then report the `Pain` and average `Pulse` in decreasing order of average `Pulse`.

- **4f (4 points).** For horses with `Pain` at least 3 and irregular `Feces` (1 means regular, 2 through 4 mean irregular), produce a plot of `Temp_Rect` versus `Pulse`, then report a cross-table of `Ab_Dist` versus `Outcome`. (By this, we mean a table that counts the number of horses with each combination of `Ab_Dist` and `Outcome` levels. To check your work, there should be 11 horses with `Ab_Dist` equal to 1 and `Outcome` equal to 1.) Hint: use the "tee" operator `%T>%`.

- **Challenge (10 points).** A test data set, on 68 additional horses, is available at https://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.test. On the original data set from Q4a, call this the training data set, fit a suitable model to predict `Ab_Dist` (abdominal distention) from the other variables. You should treat `Ab_Dist` as a categorical variable, not as a continuous variable. Use your fitted model to make predictions on the test data set, and report the percentage of correctly classified horses (the percentage of horses in the test set for which your model predicted the correct level of `Ab_Dist`). Compare this to the percentage of correctly classified horses from the "naive" model, which simply predicts the most common level of `Ab_Dist` that was observed in the training set. Does your model do significantly better? Note: missingness here is a challenge to deal with. If you were to simply remove all horses in the training set for which that had any missing data, then you would be left with 6 total (recall, you started with 300) ... so this is not a good idea. A better idea is to replace all missing values by the mean of that variable (this is a simple form of *imputation*). 