```{r echo = FALSE, results = 'hide', message = FALSE}
library(tidyverse)
library(mosaic)
library(ggthemes)
theme_set(theme_solarized())                 # change theme for ggplot2
trellis.par.set(theme=theme.mosaic()) # change default color scheme for lattice

knitr::opts_chunk$set(
    message = FALSE,
    warning = FALSE,
    tidy    = FALSE)
```

## Loading packages and data

```{r}
library(tidyverse)
library(mosaic)
library(Stat2Data)
```

```{r}
data(Leukemia)
?Leukemia
Leukemia
```

## Splitting into Train and Test

We are trying to select a model to predict `Resp` (will a Leukemia patient respond to treatment? (`1` = yes, `0` = no)) from a set of diagnostic and vitals results.

Before we do any modeling, let's set aside a test set.

```{r}
set.seed(29747) # Don't forget to do this!
TestSet <- Leukemia %>%
  slice_sample(prop = 0.20)
TrainSet <- Leukemia %>%
  setdiff(TestSet)
```

## Potential Predictors

We'll restrict our models to the following predictors:

* `Age` (age in years)
* `Temp` (highest temperature in 10 * degrees F)
* `Smear` (differential percentage of "blasts")
* `Infil` (percent of marrow leukemia infiltrate)
* `Index` (percentage labeling index of bone marrow leukemia cells)
* `Blasts` (absolute number of blasts in thousands)

**Note**: it's not a good idea to jump into modeling in a domain that you don't have any background in, especially when some of the variables are quite technical in nature.  

I chose this example because it has a little bit of "messiness", but not too much; that said, you wouldn't want to do this for a project unless you know about Leukemia, or at least about these metrics.


## Checking for potential transformations

We'll check individual (quantitative) variables' distributions first to see whether any obviously call for a transformation

### Age

```{r}
gf_boxplot(~Age, data = TrainSet) 
gf_histogram(~Age, data = TrainSet, binwidth = 10)
```

The `Age` variable is fairly symmetric.

### Temperature

```{r}
gf_boxplot(~Temp, data = TrainSet)                
gf_histogram(~Temp, data = TrainSet, binwidth = 10)
```

The `Temp` variable is kind of right-skewed. However, due to the nature of the scale (no natural "zero point") the usual transformations that we might use with a right-skewed variable don't make a lot of sense: A log isn't sensible, because ratios of degrees Fahrenheit are meaningless.

On balance I think it's probably better to leave `Temp` as is, because any transformation would be somewhat arbitrary, and sacrifice interpretability.


### Smear

```{r}
gf_boxplot(~Smear, data = TrainSet)               
gf_histogram(~Smear, data = TrainSet, binwidth = 5)
```

This one isn't too bad.


### Infil

```{r}
gf_boxplot(~Infil, data = TrainSet)
gf_histogram(~Infil, data = TrainSet, binwidth = 5)
```

Also pretty symmetric.

### Index

```{r}
gf_boxplot(~Index, data = TrainSet)               
gf_histogram(~Index, data = TrainSet, binwidth = 1)
```

That also looks pretty symmetric.


### Blasts

```{r}
gf_boxplot(~Blasts, data = TrainSet)              
gf_histogram(~Blasts, data = TrainSet)           
```


That's a lot of skew. It looks structural, though there may be one genuine outlier.

Let's remind ourselves what this variable is measuring:

```{r}
# ?Leukemia
Leukemia
```

It's a count, so ratios make sense. We could consider a log transformation, but that will produce infinite values if there are any zeroes. Let's check (I'm going to use the full dataset here, because if we do a log transform and then there are zeroes in the test set, it will break).

```{r}
fivenum(~Blasts, data = Leukemia)
tally(~(Blasts == 0), data = Leukemia)
```

There's one zero in the data, so we need to account for this if we take a log. 

One thing we could do is add a small count to each value before taking the log. The variable is measured in thousands, and appears to be recorded to one decimal place, which means the counts are rounded to the nearest hundred. So let's add 0.05, which is the equivalent of 50 blasts, which is the amount of potential "rounding error" inherent in each measurement.

The base of the log doesn't affect anything except the scaling of the result, so I'll use base 2 so that one unit represents a doubling of the count, which seems like a natural thing to think about with cells that are dividing.

```{r}
gf_boxplot(~log2(Blasts + 0.05), data = TrainSet)
gf_histogram(~log2(Blasts + 0.05), data = TrainSet, binwidth = 0.5)
```


That looks a lot better. That high value probably isn't an outlier after all, but rather just a representative of a long right tail in the distribution.


## Our pool of variables

All of the candidate predictors had pretty symmetric distributions, with the exception of `Temp`, which we elected not to transform because there wasn't a natural transformation to do (and the skew wasn't all that bad to begin with), and `Blasts`, which we decided to transform via base 2 `log`, after adding a count equal to the amount of potential rounding error in the original variable.

* Age --> Age (no transformation)
* Temp --> Temp (no transformation)
* Smear --> Smear (no transformation)
* Infil --> Infil (no transformation)
* Index --> Index (no transformation)
* Blasts --> log2(Blasts + 0.05)


## Checking for multicollinearity

Next, let's check whether any of the variables are redundant with the others.

```{r, message = FALSE}
library(car)
## Since vif() expects a model, we'll fit a model with all predictors,
## only for the purposes of checking multicollinearity
fullModel <- 
  glm(Resp ~ Age + Temp + Smear + Infil + Index + log2(Blasts + 0.05),
     family = "binomial", data = TrainSet)
vif(fullModel)
```

Smear and Infil are fairly predictable from the others, but neither exceeds our VIF > 5 guideline. We won't preemptively remove anything, but it's worth keeping in mind that we have some overlap in information.


## Choosing combinations of predictors to assess.

Age and Temp are easy variables to collect, so let's start with a baseline model that uses those and none of the biomarker variables, and then consider which of the biomarkers are the most important to include.

As a quick first check, let's see if any of them are useful:

```{r}
baseModel <- glm(
  formula = Resp ~ Age + Temp, 
  family  = "binomial", 
  data    = TrainSet)
fullModel <- glm(
  formula = Resp ~ Age + Temp + Smear + Index + Infil + log2(Blasts + 0.05), 
  family  = "binomial", 
  data    = TrainSet)
```

```{r}
anova(baseModel, fullModel, test = "LRT")
```

There's fairly strong evidence that at least some of those biomarkers are helpful, after controlling for `Age` and `Temp`.


## Assessing models using AIC

Squared deviations aren't something we use with a binary response, so $R^2$ isn't a relevant measure. 

Instead, we can look at AIC, which stands for "Akaike Information Criterion", which is analogous to Mallow's $C_p$, and is a measure of how well a model fits a dataset which is based on the log likelihood, and includes a penalty for complexity. Like Mallow's Cp, smaller values are better.

Since there are four markers in question, there are sixteen (= 2^4) possible combinations of these markers we could include in our model. We could fit all of these and examine their AIC, but to cut down on the number we consider, let's first see if there is one of these markers that provides most of the information we need on its own.

```{r}
indexModel <- glm(
  formula = Resp ~ Age + Temp + Index,
  family  = "binomial", 
  data    = TrainSet)
smearModel <- glm(
  formula = Resp ~ Age + Temp + Smear,
  family  = "binomial", data = TrainSet)
infilModel <- glm(
  formula = Resp ~ Age + Temp + Infil,
  family  = "binomial", data = TrainSet)
blastsModel <- glm(
  formula = Resp ~ Age + Temp + log2(Blasts + 0.05),
  family  = "binomial", 
  data    = TrainSet)
```

```{r}
AIC(smearModel)
AIC(blastsModel)
AIC(infilModel)
AIC(indexModel)
AIC(fullModel)
```

Including `Index` marker and no other biomarkers actually gives us an AIC which is better than the one we get using all of the biomarkers.

Because Index yields such a good AIC all by itself, we likely want it included. 

Next, let's fit the various models with Index and one other marker, and see if any of these can improve on the AIC we get from the index-only model.

```{r}
indexSmearModel <- glm(
  formula = Resp ~ Age + Temp + Index + Smear,
  family  = "binomial", 
  data    = TrainSet)
indexInfilModel <- glm(
  formula = Resp ~ Age + Temp + Index + Infil,
  family  = "binomial", 
  data    = TrainSet)
indexBlastModel <- glm(
  formula = Resp ~ Age + Temp + Index + log2(Blasts + 0.05),
  family = "binomial", 
  data   = TrainSet)
```

```{r}
AIC(indexSmearModel)
AIC(indexInfilModel)
AIC(indexBlastModel)
AIC(indexModel)
AIC(fullModel)
```

Looks like `Index` and `Infil` together give us a potentially better model (at least, when fitting on this dataset) than including all four markers. So do `Infil` and `Smear` together. What about all three?

```{r}
indexSmearInfilModel <- glm(
  formula = Resp ~ Age + Temp + Index + Smear + Infil,
  family  = "binomial", 
  data    = TrainSet)
```

```{r}
AIC(indexModel)
AIC(indexSmearModel)
AIC(indexInfilModel)
AIC(indexBlastModel)
AIC(indexSmearInfilModel)
AIC(fullModel)
```

Including all three yields a higher AIC than including index and either one of the others. This is consistent with our previous observation that `Smear` and `Infil` are somewhat multicollinear with other predictors, so it may not be worthwhile to include both together. However, all of these models give lower AICs than the full model.

Let's check the residual plots for these models.

```{r}
library(arm)
binnedplot(
  x = fitted(indexModel),
  y = residuals(indexModel, type = "pearson"))
binnedplot(
  x = fitted(indexSmearModel),
  y = residuals(indexSmearModel, type = "pearson"))
binnedplot(
  x = fitted(indexInfilModel),
  y = residuals(indexInfilModel, type = "pearson"))
binnedplot(
  x = fitted(indexBlastModel),
  y = residuals(indexBlastModel, type = "pearson"))
binnedplot(
  x = fitted(indexSmearInfilModel),
  y = residuals(indexSmearInfilModel, type = "pearson"))
binnedplot(
  x = fitted(fullModel),
  y = residuals(fullModel, type = "pearson"))
```

None of them show particular evidence of curvature.


## Checking for possible interactions

One question we might have is whether the relationship between any of the biological markers and the response to treatment is different for patients of different ages.  That is, should we interpret the markers differently depending on the age of the person? 

This is equivalent to asking whether we want to include an interaction term between Age and any of the other variables.

```{r}
ageInteractionsModel <- glm(
  formula = Resp ~ Temp + (Index + Infil + Smear + Blasts) * Age,
  family  = "binomial", 
  data    = TrainSet)
anova(fullModel, ageInteractionsModel, test = "LRT")
AIC(ageInteractionsModel)
```

There's not much evidence that including interaction terms with `Age` helps the fit, and the AIC is worse.


## Cross-Validation

We've settled on a pool of six models. All six include `Temp`, `Age` and `Index`. Then we also include

1. Nothing else
2. `Infil`
3. `Smear`
4. `Blasts`
5. `Infil` and `Smear`
6. `Infil`, `Smear` and `Blasts`

Let's do cross-validation with these six models.

```{r}
# source("http://colindawson.net/stat213/code/helper_functions.R")
source("~/www/colinreimerdawson/stat213/code/helper_functions.R")
```

Some error metrics that make sense for logistic regression:

**"Hard" misclassification rate**: Count prediction as an "error" if fitted probability is on the "wrong side" of 0.5 and compute proportion of validation cases that are "misclassified".

* Can use `hard_classification_error()` from `helper_functions.R`

**False positive rate:**: What proportion of the true 0s are misclassified as 1s (this ignores the true 1s)

* Use `false_positive_rate()`

**False negative rate:**: What proportion of real 1s are misclassified as 0s (ignores true 0s)

* Use `false_negative_rate()`

**"Soft" misclassification rate:**: Count each point as a "partial error" according to the predicted probability assigned to the opposite category.

* Rewards confident correct predictions and penalizes confident wrong predictions, with "uncertain" predictions (near 0.5) in between.
* Computes average error per case in validation set
* Use `soft_classification_error()`

**Normalized deviance**: Residual deviance (on validation set), divided by sample size (i.e., "per point" deviance)

* Use `normalized_deviance()`

**Concordance, discordance, and Somers' D:**

* These measures consider all pairs of cases with distinct response values (that is, one 0 and one 1)
* Count what proportion of the time the 1 is assigned a higher predicted prob. than the zero (good).  This is the **"concordance"**.
* Count what proportion of the time the opposite occurs. This is the **"discordance"** (bad)
* Note that these may not add to 1, since if the predicted probability is the same, then it counts toward neither measure
* **Somers' D** := (Concordance - Discordance) / TotalPairs

We'll make a list of our six finalists, plus the "base" model (for comparison).
```{r}
finalists <- list(
  IndexOnly       = indexModel, 
  IndexInfil      = indexInfilModel,
  IndexSmear      = indexSmearModel,
  IndexBlasts     = indexBlastModel,
  IndexSmearInfil = indexSmearInfilModel,
  AllFour         = fullModel,
  None            = baseModel)
formulaList = sapply(finalists, formula)
```

Now we'll do 10-fold cross-validation, with 10 random splits for each (to reduce the impact of the choice of folds):

```{r}
set.seed(8675309)
normalized_deviances <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = normalized_deviance, 
  iterations     = 10,
  family         = "binomial")
normalized_deviances
```

All six are relatively close compared to the gap between them and the "base" model, but `IndexSmear` has the lowest normalized deviance.

```{r}
set.seed(8675309)
## hard misclassification rate (lower is better)
hard_misclassification_rates <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = hard_classification_error, 
  iterations     = 10,
  family         = "binomial")  
hard_misclassification_rates
```

The model with only `Index` and `Smear` again has the best result.

Let's examine the false-positive and false-negative rates separately, since if we were using this model to decide whether treatment was worthwhile, deciding to treat and then getting no response presumably has a different "cost" than deciding not to treat if there would have been a response.

```{r}
set.seed(8675309)
false_positive_rates <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = false_positive_rate, 
  iterations     = 10,
  family         = "binomial")    
false_positive_rates
```

This time, the model with Index and Blasts gives the best result.

What about false negatives?

```{r}
set.seed(8675309)
false_negative_rates <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = false_negative_rate, 
  iterations     = 10,
  family         = "binomial")
false_negative_rates
```

`Index` and `Smear` yields the lowest false negative rate, however.

How about "soft misclassification rate"?

```{r}
set.seed(8675309)
soft_misclassification_rates <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = soft_classification_error,
  iterations     = 10,
  family         = "binomial")
soft_misclassification_rates
```

`Index` and `Infil` gives a slightly lower soft misclassification rate than `Index` and `Smear`, but it's really close.

```{r}
set.seed(8675309)
## Higher is better
concordances <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = concordance,
  iterations     = 10,
  family         = "binomial")
concordances
```

The model with `Index` and `Infil` performs the best on this measure.

```{r}
set.seed(8675309)
## Lower is better
discordances <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = discordance,
  iterations     = 10,
  family         = "binomial")
discordances
```

and on this one.

```{r}
set.seed(8675309)
## Higher is better
somers_Ds <- sapply(
  X              = formulaList,
  FUN            = do_k_fold_cv, 
  dataset        = TrainSet, 
  K              = 10,
  error_function = somers_d,
  iterations     = 10,
  family         = "binomial")
somers_Ds
```

Since Somer's d is calculated by combining concordance and discordance, it makes sense that the model that performed best on both of these separately performs the best here.

## Conclusion

Other than false positive rate, which was minimized by using `Blasts`, two models consistently perform the best: `Index` and `Infil`, and `Index` and `Smear`.

Which of these we would prefer depends on context; e.g., how invasive and expensive is the procedure to collect the `Smear` variable compared to the `Infil` variable, how heavily we want to weigh false negatives vs being able to "rank" response probabilities, etc.


## Evaluating performance on the test set.

Whichever model we decide to go with, we should calculate our various performance measures on the held out test set.

```{r}
metric_list <- list(
  "Hard Misclassification" = hard_classification_error,
  "Soft Misclassification" = soft_classification_error,
  "Normalized Deviance"    = normalized_deviance,
  "False Positive Rate"    = false_positive_rate,
  "False Negative Rate"    = false_negative_rate,  
  "Concordance"            = concordance,
  "Discordance"            = discordance,
  "Somer's D"              = somers_d)

test_metrics <- sapply(
  X            = metric_list,
  FUN          = held_out_error,
  test_dataset = TestSet,
  fitted_model = indexInfilModel)
test_metrics
```


## Fit on the full dataset and interpret

```{r}
finalModel <- 
  glm(
    formula = Resp ~ Age + Temp + Index + Infil, 
    family  = "binomial", 
    data    = Leukemia)
finalModel %>% 
  coef() %>% 
  round(3)
```

### Prediction equation:

$$ \widehat{LogOdds(Resp = 1)} = 95.568 - 0.060\cdot Age - 0.099 \cdot Temp + 0.407 \cdot Index + 0.034 \cdot Infil$$


### Odds and odds ratios:

```{r}
finalModel %>% 
  coef() %>% 
  exp() %>% 
  round(3)
```

```{r}
Leukemia
```


It seems that each year of Age is associated with about a 5.8% decrease in the odds of responding to treatment, controlling for temperature and biomarkers. 

Patients with higher temperature readings have lower odds of responding, about 9.5% lower for each tenth of a degree Fahrenheit. 

Higher `Index` scores are associated with a higher probability of responding to treatment, after controlling for age, temperature, and percentage infiltrate, with the odds of response increasing by about 50% with each 1 point increase in Index relative to what is expected based on the other predictors.

Similarly, higher `Infil` scores are associated with a higher probability of responding to treatment, after controlling for the other predictors, with the odds of response increasing by about 3.5% for each 1 point difference.


## Confidence intervals

Let's look at the predicted recovery probabilities for a few example patients.

```{r}
f_hat <- makeFun.logistic(finalModel)
```

To see what sensible values to plug in for `Index` and `Infil` are, let's look
at a scatterplot and five number summaries for these markers.

```{r}
fivenum(~Index, data = Leukemia)
fivenum(~Infil, data = Leukemia)
```

```{r}
gf_point(Index ~ Infil, data = Leukemia)
```


**Patient 1a: 20 years old with relatively high index and infil scores, and a normal body temperature** (favorable scenario)

```{r}
f_hat(Age = 20, Temp = 986, Index = 13.5, Infil = 72.5, interval = "confidence")
```

We estimate that a patients with high index and infil scores and a normal body temp have at least an 83% chance of responding to treatment (the lower endpoint of the confidence interval)

**Patient 1b: same as Patient 1a but with a fever of 101F**

```{r}
f_hat(Age = 20, Temp = 1010, Index = 13.5, Infil = 72.5, interval = "confidence")
```

If a patient like this has a fever, the probability of responding might still be as high as 97 or 98%, 
but it might be down to something like 51%.

**Patient 2a: 50 years old, moderate index and infil scores (near the median), normal body temp**

```{r}
f_hat(Age = 50, Temp = 986, Index = 9, Infil = 61, interval = "confidence")
```

Patients like this might have anywhere between a 34% and 75% chance of responding; our best guess is that their response chance is around 56%.

**Patient 2b: Same as 2a but with a fever of 101F**

```{r}
f_hat(Age = 50, Temp = 1010, Index = 9, Infil = 61, interval = "confidence")
```

The presence of a fever of 101F worsens the outlook for this patient considerably. The probability of responding may be as low as 2%, and our best guess is that it's around 10%.

**Patient 3: 75 years old, fever of 101F, low index and infil (unfavorable scenario)**

```{r}
f_hat(Age = 75, Temp = 1010, Index = 6, Infil = 37.5, interval = "confidence")
```

Even if our model is pessimistic, we doubt that the probability of responding to treatment 
is higher than 10% or so for a patient like this, and more likely it's less than 1%.