```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("~/Github/bigdata/lectures/l7/")
library(tidyverse)
library(glmnet)
```

## Today's Class {.smaller}

1. Review
    - LASSO
    - Deviance vs MSE
    - Cross-Validation
    - R help docs -- An explainer
2. KNN
3. Binomial Classificiation: Probabilities to Predictions
4. Misclassification: Sensitivity, Specificity, ROC Curve
5. Multinomial Regression

# Quick Review

## LASSO

 The lasso penalty is the $l_1$ norm of our parameters: $pen(\beta) = \sum |\beta_j|$, which penalizes larger coefficients and non-zero coefficients. So our estimates are:
<span style="color:maroon">
$$\hat{\beta}_\lambda = \underset{\beta}{\text{argmin}}\left( Dev(\beta)+\lambda \sum_{j=1}^p |\beta_j|\right)$$</span> 

For a given model of the error terms, $Dev(\beta)$ is held fixed, so we can make comparisons within the logistic or linear regressions, but not between them easily. 

When we switch logistic $\rightarrow$ linear, we change $Dev$. We can still make comparisons betwen these models but we need to ensure we are comparing the same thing.

## LASSO -- MSE vs Deviance

For the logistic LASSO, our deviance is approximately:
$$Dev_{logit}(\beta) \propto \sum_{i=1}^n [ log(1+exp(x_i'\beta))-y_ix_i'\beta]$$

While for our linear LASSO, our deviance is approximately:
$$Dev_{linear}(\beta) \approx  {1 \over \sigma^2} \sum_{i=1}^n (y-x'\beta)^2 \propto \widehat{MSE}$$

This difference derives from the different *likelihoods* each model uses to model their different ideas about the error distribution. 

<span style="color:maroon">But there is only one true error distribution.</span> 

## Model Choice
When deciding which model to use, we used the out-of-sample deviance to select a value of $\lambda$ for each LASSO type. But we still need to choose between these two very different models. 

<span style="color:darkgreen">That choice is still driven by out of sample prediction errors</span>. <span style="color:maroon">We need to compare the prediction errors using the same yardstick.</span> 

So we could compare all the prediction errors using the binomial deviance, or using the MSE, or using some other loss function. 

## Model Choice

Critically (and we will see this again), cross-validation gives us a good understanding of our out-of-sample error. 

We can answer questions like <span style="color:lightblue"> 
 "When I build a model in this way, what do its out of sample errors look like?"</span>
 
Which means we can also answer questions like <span style="color:maroon">"Which model building process has better out-of-sample predictions?" </span>  

All of this is because of the very trustworthy error distributions we get from Cross-validation. 

## R help documents. 

I've realized I should give you a quick overview of how the R help documents are structured. 

There are three main questions I turn to the help documents with.

1. Where in the output is $x$? 
2. What do I put in the input to do $y$?
3. How do I do $z$?

By far, the help documents are most useful for (1) and (2). 

```{r, eval=F}
?cv.glmnet
```

## Basic Structure of Help Docs {.smaller}

> 1. **Description**: Basic explanation of the function in question. 
> 2. **Usage**: This shows the function (and close relatives), as well as most of the arguments to the function, and their defaults. 
> 3. **Arguments**: This lists out every argument you can enter, describing what it must be, and what it controls. 
> 4. **Details**: In-depth description of the function, often including further detail about inputs, sometimes math.
> 5. **Value**: In depth description of the output of the function. 
> 6. **Authors/Refs**: People and places with even more details
> 7. **See Also**: other closely related functions. Many model tools will mention "`plot`, `predict`, and `coef`" methods for `cv.glmnet`. That implies functions like `coef.cv.glmnet` and `plot.cv.glmnet` you can look up.
> 8. **Examples**: example code running the function on basic data. 

# Classification

## Basic Setting
Just as in our basic prediction problems, we have data with $n$ observations $(\textbf{x}_i,y_i)$ of something. 

But now $y_i$ is qualitative rather than quantitative. Membership in some category $\{1,...,M\}$

The basic problem then is the following: <span style="color:maroon"> 
Given new observation covariates $\textbf{x}_i^{new}$, what is the class label $y_i^{new}$? </span>

The quality of any classifier can be determined by its misclassification risk:
$$P[\hat{y}_i^{new} \neq y_i^{new}]$$

## Motivation

How does this differ from basic logistic question of predicting $P[y_i^{new}=1]$?

1. We may have many categories, not just two.
2. We may have different discrete actions we take depending on classification. 
    - In some domains, as $P[y_i=1]$ changes, our actions change smoothly. 
        - e.g. as $P[$stock x goes up$]$ we buy more of it.
    - In other domains, as $P[y_i=1]$ changes, our actions change suddenly. 
        - e.g. as $P[$get into college i$|$I apply$]$ goes up, we switch from "don't apply" to "apply", there is mostly no "apply a little bit more"

## Motivation

This leads to situations where we want to categorize, as it affects our decisions. 

We face decisions not on a spectrum, so the most useful interpretation of our predictions may not be on that spectrum. 

$\implies$ classification. 

## Optimal Classifier

Presuming you have no preference between different types of misclassification (**LOSS FUNCTION CLAIM**), there is an optimal classifier, known as the <span style="color:maroon">Bayes Classifier</span>. 

$$ \hat{y}_i^{new} = \underset{j \in \{1,...,M\}}{\text{argmax}} P[y_i^{new} = j |\textbf{x}_i^{new}]$$

Find the prediction which is *most* likely (not necessarily all that likely -- e.g. $P[\hat y = y] < 0.5$ is common). This will minimize the misclassification risk.

## Bayes Classifier

Unfortunately, we don't know $P[y = j | \textbf{x}_i^{new}]$. So the Bayes classifier is in some ways an unattainable standard. 

But we can estimate it!

## Estimating P

There are many tools for estimating $P[y|x]$ given the training data. 

- We can use parametric tools.
    - Assume $P[y|x]$ is some function of unknown parameters $\beta$, estimate those, and make predictions.
        - Sound familiar? <span style="color:lightblue">Logistic Regression</span> does this.
- We can use non-parametric tools. 
    - Estimate $P[y|x]$ directly without any parameters. 
        - <span style="color:maroon">K Nearest Neighbors</span> (KNN)

# KNN

## KNN Basics

**Basic Idea**: Estimate $P[y|x]$ locally using the labels of similar observations in the training data. 

KNN: What is the most common class near $x^{new}$? 

1. Take the $K$ nearest neighbors $x_{i,1},...,x_{i,K}$ of $x^{new}$ in the training data
    - Nearness is (usually) Euclidean distance: $\sqrt{\sum_{j=1}^p (x^{new}_j-x_{i,k,j})^2}$
2. Estimate $P[y=j|x] = \sum_{i=1}^K 1(y_i=j)$
3. Select the class $j$ with the highest probability.

## KNN -- Details

This (again) is sensitive to the scale of each covariate $x$. So we will rescale them all by standard deviations. 

This will be sensitive to $K$, and we need to pick it. 

- $K=n$ -- we just take the mean across the entire training data.
- $K=1$ -- Whatever observation happens to be closest will be our prediction.

Cross-validation here will help. 

## KNN Example Data

```{r, echo=F}
n = 50
set.seed(100)
set.seed(101)
x1 = runif(n)
x2 = runif(n)
prob  = ifelse(x1 < 0.5 & x1 > 0.25 & x2 > 0.25 & x2<0.75,0.8,0.3)
newx = data.frame(y="New Observation",x1=0.7,x2=0.6)
y  = as.factor(rbinom(n,1,prob))
levels(y) = c("1","2","New Observation")
df = data.frame(y=y,x1=x1,x2=x2)
df = rbind(df,newx)
```

```{r,echo=F}
ggplot(df,aes(x=x1,y=x2,col=y)) +
  geom_point()+
  theme(aspect.ratio = 1)
```

## KNN Example $K=3$

```{r,echo=F}
dists = sqrt((x1-0.7)^2+(x2-0.6)^2)
k3 = sort(dists)[3]
gg_circle <- function(r, xc, yc, color="black", fill=NA, ...) {
    x <- xc + r*cos(seq(0, pi, length.out=100))
    ymax <- yc + r*sin(seq(0, pi, length.out=100))
    ymin <- yc + r*sin(seq(0, -pi, length.out=100))
    annotate("ribbon", x=x, ymin=ymin, ymax=ymax, color=color, fill=fill, ...)
}
ggplot(df,aes(x=x1,y=x2,col=y))+
  gg_circle(k3,0.7,0.6,fill="salmon",alpha=0.3,lwd=0)+
  geom_point()+
  theme(aspect.ratio = 1)
```

## KNN Example $K=7$

```{r,echo=F}
k7 = sort(dists)[7]
ggplot(df,aes(x=x1,y=x2,col=y))+
  gg_circle(k7,0.7,0.6,fill="green",alpha=0.3,lwd=0)+
  geom_point()+
  theme(aspect.ratio = 1)
```
The relative 'vote counts' are a very crude estimate of probability.

## KNN Example $K=1$


```{r, echo=F}
k1 = sort(dists)[1]
ggplot(df,aes(x=x1,y=x2,col=y))+
  gg_circle(k1,0.7,0.6,fill="salmon",alpha=0.3,lwd=0)+
  geom_point()+
  theme(aspect.ratio = 1)
```

## KNN Example

```{r,echo=F}
grid.fineness = 81
sequence = seq(0,1,length.out=grid.fineness)
grid = expand.grid(sequence,sequence)
knn_pred = function(x,k=5) {
  dists = sqrt((x1-x[1])^2+(x2-x[2])^2) #Find all distances to current obs
  bound = sort(dists)[k]                #Find kth smallest distance
  indices = which(dists <= bound)       #Find which obs have dists 1:k
  outcomes = as.integer(y[indices])     #Find corresponding outcomes y
  round(mean(outcomes)) #Taking advantage of 2 outcomes. 
}
colnames(grid) = c("x1","x2")
yhat = as.factor(apply(grid,1,knn_pred))
levels(yhat) = c(levels(yhat),"New Observation")
df_grid = data.frame(x1=grid$x1,x2=grid$x2,y=yhat)
ggplot(df,aes(x=x1,y=x2,col=y))+
  geom_point(data=df_grid,aes(x=x1,y=x2,col=y),alpha=0.5,size=0.2)+
  geom_point()+
  theme(aspect.ratio=1)
```


## KNN Example: More Data

```{r,echo=F}
n = 2000
x1 = runif(n)
x2 = runif(n)
prob  = ifelse(x1 < 0.5 & x1 > 0.25 & x2 > 0.25 & x2<0.75,0.8,0.3)
newx = data.frame(y="New Observation",x1=0.7,x2=0.6)
y  = as.factor(rbinom(n,1,prob))
levels(y) = c("1","2","New Observation")
df = data.frame(y=y,x1=x1,x2=x2)
df = rbind(df,newx)
dists = sqrt((x1-0.7)^2+(x2-0.6)^2)
ggplot(df,aes(x=x1,y=x2,col=y))+geom_point(alpha=0.5)+theme(aspect.ratio=1)
```



## KNN Example $K=3$
```{r,echo=F}
grid.fineness = 101
sequence = seq(0,1,length.out=grid.fineness)
grid = expand.grid(sequence,sequence)
colnames(grid) = c("x1","x2")
yhat = as.factor(apply(grid,1,knn_pred,k=3))
levels(yhat) = c(levels(yhat),"New Observation")
df_grid = data.frame(x1=grid$x1,x2=grid$x2,y=yhat)
ggplot(df_grid,aes(x=x1,y=x2,col=y))+
  geom_point(size=0.2)+
  theme(aspect.ratio=1)
```


## KNN Example $K=11$

```{r, echo=F}
yhat = as.factor(apply(grid,1,knn_pred,k=11))
levels(yhat) = c(levels(yhat),"New Observation")
df_grid = data.frame(x1=grid$x1,x2=grid$x2,y=yhat)
ggplot(df_grid,aes(x=x1,y=x2,col=y))+
  geom_point(size=0.2)+
  theme(aspect.ratio=1)
```






## KNN Example $K=31$
```{r,echo=F}
yhat = as.factor(apply(grid,1,knn_pred,k=31))
levels(yhat) = c(levels(yhat),"New Observation")
df_grid = data.frame(x1=grid$x1,x2=grid$x2,y=yhat)
ggplot(df_grid,aes(x=x1,y=x2,col=y))+
  geom_point(size=0.2)+
  theme(aspect.ratio=1)
```

## KNN -- Example Base Truth
The optimal prediction scheme

```{r,echo=F}
opt_pred = function(x){
  out = 1
  if(x[1] < 0.5 & x[1]>0.25 & x[2] < 0.75 & x[2] > 0.25)
    out = 2
  out
}
opt_preds = apply(grid,1,opt_pred)
df_grid_opt = df_grid
df_grid_opt$y = as.factor(opt_preds)
ggplot(df_grid_opt,aes(x=x1,y=x2,col=y))+
  geom_point(size=0.2)+
  theme(aspect.ratio=1)
```

## KNN

Higher $K$ leads to higher *in-sample* training error. (Proportion incorrect).

Lower $K$ leads to higher flexibility $\implies$ overfitting and poor OOS misclassification.


## KNN {.smaller}

- Pros:
    - Straightforward
    - Easily adjust to multiple categories.
    - Outperform other classifiers in settings where local details matter
- Cons:
    - Computing neighbors (just finding them) can be costly with big $n$ or $p$
    - No variable selection built in
    - Choosing $k$ is tricky -- Cross-validation can work -- but it is unstable.
        - Remember, LASSO benefitted from *stability* of estimates
    - Classification is very sensitive to $k$
    - Only other output is rough local probabilities.
        - Without good probabilities, assessing uncertainty is hard.
        
## Binary Classification

Many problems can be reduced to binary classification (as above). 

KNNs are a useful non-parametric classification tool.

*Logits* are a useful parametric classification tool.  
- Remember Spam?

- Logits yield parametric decision boundaries. Easy to interpret.
- Logits are global methods. Use *all* the training data to inform predictions.
    - The probability estimates are more useful and more stable.
- Logits can do variable selection.

## Credit Classification Example

German loan/default data. "Predict performance of new loans" -- want to predict default probability. 

Going to try to use borrower and loan characteristics. 

Messy Data.

## Beware
This is not 'randomly sampled' data.

```{r, echo=F}
## read data and create some `interesting' variables
credit = read.csv("~/Github/bigdata/lectures/l7/credit/credit.csv")
## re-level the credit history and checking account status
credit$history = factor(credit$history, levels=c("A30","A31","A32","A33","A34"))
levels(credit$history) = c("good","good","poor","poor","terrible")
## a few others

credit$foreign = factor(credit$foreign, levels=c("A201","A202"), labels=c("foreign","german"))
credit$rent = factor(credit$housing=="A151")
credit$purpose = factor(credit$purpose, levels=c("A40","A41","A42","A43","A44","A45","A46","A47","A48","A49","A410"))
levels(credit$purpose) = c("newcar","usedcar",rep("goods/repair",4),"edu",NA,"edu","biz","biz")
```

```{r, echo=F}
plot(factor(Default) ~ history, data=credit, col=c(8,2), ylab="Default") ## surprise!
```

## Beware
Consider your data sources carefully. 

```{r, echo=F}

plot(factor(Default) ~ purpose, data=credit, col=c(8,2), ylab="")
```


## German Credit LASSO
```{r, echo=F}
xnaref = function(x){
	if(is.factor(x))
		if(!is.na(levels(x)[1]))
			x = factor(x,levels=c(NA,levels(x)),exclude=NULL)
	x
}
naref = function(DF){
	if(is.null(dim(DF))) return(xnaref(DF))
	if(!is.data.frame(DF)) 
		stop("You need to give me a data.frame or a factor")
	DF = lapply(DF, xnaref)
	return(as.data.frame(DF))
}

credit = naref(credit)
credx = sparse.model.matrix( ~ (duration + amount +
    installment + age + history +
    purpose + foreign + rent)^2, data=credit)
default = credit$Default
```
```{r}
credscore = cv.glmnet(credx, default, family="binomial")
plot(credscore)
```

## German Credit LASSO

```{r}
sum(coef(credscore, s="lambda.1se")!=0) # 1se
sum(coef(credscore, s="lambda.min")!=0) # min
```

## Decision Making

There are two ways to be wrong in this binary problem.

- <span style="color:maroon">False Positive</span> predict $\hat y=1$ when $y=0$. 
    - Classify as defaulters when they are not
- <span style="color:maroon">False Negatives</span> predict $\hat y=0$ when $y=1$.
    - Classify as non-defaulters when they are.
    
Both mistakes are bad, but sometimes one is worse than the other. Logistic regression gives us an estimate of $P[y=1|x]$. And the Bayes decision rule classifies purely based on probabilities. <span style="color:lightblue">When $P[y=1|x] > 0.5$, classify as 1.</span> 
 
But instead of minimizing misclassification risk, we want to minimize our **loss**. 

## Using probabilities for Decisions

To make optimal decisions you need to account for probabilities and costs.

If for each loan, you make \$0.25 when repaid, and lose \$1 when not, then we only *expect* to make 'profits' when $P[y=1] < 0.2$. 

So we may want our classifier to use a different threshold. 

## False Positive Rate

Much like in lecture 1, we we can think about our *rate* of false positives. 

- False Positive Rate  = #Misclassified as 1 / #classified as 1
- False Negative Rate =  #Misclassified as 0/ #Classified as 0

## Sensitivity and Specificity

But we may also want to think about sensitivity and specificity.

- Sensitivity: proportion of true $y=1$ classified as such.
- Specificity: proportion of true $y=0$ classified as such.

A rule is sensitive if it mostly gets the 1s right. A rule is specific if it mostly gets the 0s right.

## ROC Curve

We can plot the ROC curve for different choices of threshold. 

```{r, echo=F}
pred = predict(credscore, credx, type="response")
pred = drop(pred) 

roc = function(p,y, ...){
  y = factor(y)
  n = length(p)
  p = as.vector(p)
  Q = p > matrix(rep(seq(0,1,length=100),n),ncol=100,byrow=TRUE)
  specificity = colMeans(!Q[y==levels(y)[1],])
  sensitivity = colMeans(Q[y==levels(y)[2],])
  plot(1-specificity, sensitivity, type="l", ...)
  abline(a=0,b=1,lty=2,col=8)
}

roc(p=pred, y=default, bty="n")
```

# Wrap up

## Things to do

HW 3 is due tomorrow night.

See you Thursday

## Rehash

- Classification is a different type of problem
- KNN is a useful tool for using local data to make predictions. 
- Logits can also be used for classifications
- Prediction Errors in this setting are closely related to eachother

# Bye!








