# Logistic Regression
We will work with a dataset packaged in the library `DAAG`:

```{r }
library(DAAG)
data(frogs)
```

We will visualize this data by plotting the presence/absence of frogs geographically (`northing` is the distance north, `easting` is distance east in meters).
We first define a factor from our binary variable for coloring the points

```{r }
presAbs<-factor(frogs$pres.abs,levels=c(0,1),labels=c("Absent","Present"))
```

We can chose the `pch` based on the values of `presAbs`. `16` means to fill in the circle, while `1` leaves it empty:

```{r }
plot(northing ~ easting, data=frogs, pch=c(1,16)[presAbs],
     xlab="Meters east of reference point", ylab="Meters north")
legend("bottomleft", legend=levels(presAbs), pch=c(1,16))
```

We will also consider another dataset, `spam7` that gives classification of emails as spam or not.

```{r }
data(spam7)
head(spam7)
```

For simplicity, we are going to rename the dataset to just `spam`

```{r }
spam = spam7
```

We use the function `glm` to fit the logistic model to the frogs data. 
We will not use the geographic variables in our model, so we will first create a data.frame that excludes those.

```{r }
frogsNoGeo<-frogs[,-c(2,3)]
```

We now call `glm` using the same syntax as in `lm`. We give the response variable on the left of the `~`, and on the right we give the explanatory variables. Because we removed the geographic variables, we will use all of the remaining variables as explanatory variables, so we can simplify the formula by just using "." on the right of the `~`, like with `lm`.

```{r }
glmFrogs = glm(pres.abs ~ ., family = binomial, data = frogsNoGeo)
```

As with `lm`, the `summary` of the output gives relevant information regarding the significance of coefficients and the overall fit.

```{r }
summary(glmFrogs)
```

Like in `lm`, we can use the function `fitted` to get the fitted values of our data. These are the fitted probabilities for each observation.

```{r }
head(fitted(glmFrogs))
```

We can plot the fitted probabilities against the actual response. However, since our actual responses aren't probabilities, but 0/1 values, it makes more sense to plot this as a boxplot, so I can see the distribution of the fitted values for the 0 group and the 1 group. 

```{r }
boxplot(fitted(glmFrogs)~frogs$pres.abs, at=c(0,1))
points(x=jitter(frogs$pres.abs),fitted(glmFrogs))
```

I also overlay the individual predictions on top of the boxplot. To do this, I used `points`. But if I had just done `points(x=frogs$pres.abs, y=fitted(glmFrogs))`, then all the x values in the `0` group would be the same and they would all be overplotted. So I used the function `jitter` to add a bit of random noise to the 0/1 values; this helps spread out the x-axis values so the points don't get overlaid on each other.

We can see some of the points have high predicted probabilities to be in the opposite class from what we observed. For example, some points that had a response of 0 had fitted probabilities $>0.7$ to be a 1. Similarly some points with a response of 1 had fitted probabilities $<0.2$ to be a 1. We can isolate those points, and inspect them to try to understand why.

First we find which points these are:

```{r }
high0 <- frogs$pres.abs==0 & glmFrogs$fitted >.7
low1 <- frogs$pres.abs==1 & glmFrogs$fitted <.2
```

Once we've found them, we can plot them as red on our plot. 

```{r }
par(mfrow=c(1,2))
plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1], col=c("black","red")[factor(high0)], xlab="Meters east of reference point", ylab="Meters north",main="Points with no frogs, but high prob") 
plot(northing ~ easting, data=frogs, pch=c(1,16)[frogs$pres.abs+1], col=c("black","red")[factor(low1)], xlab="Meters east of reference point", ylab="Meters north",main="Points with frogs, but low prob") 
```

# Comparing models
We can calculate the overall residual deviance from our object using the `deviance` function:

```{r }
deviance(glmFrogs)
```

We demonstrate the difference in deviance we get from removing the variable `NoOfPools`. First we fit a model excluding this variable:

```{r }
m2 = glm(pres.abs ~ altitude + distance + NoOfSites + avrain + meanmin + meanmax, family = binomial, data = frogs) 
```

Now we calculate the deviance for this smaller model and the original (full) model:

```{r }
deviance(m2)
deviance(glmFrogs)
```

We can similarly get the deviance for a null model of no variables. This is the 

```{r }
m0 <- glm(pres.abs ~ 1, family = binomial, data = frogs)
deviance(m0) 
```

Notice that this null deviance is what is reported in the summary of our glm object:

```{r }
summary(glmFrogs)
```

This comparison of deviances can also be compared via the function `anova`

```{r }
anova(m0,glmFrogs)
```

The above did not give a p-value for the comparison, but if we specify the test we want to use, such as `test="LRT"`, you will get p-values.

```{r }
anova(m0,glmFrogs, test="LRT")
```

We can also get p-values for compare the effect of leaving out a single variable:

```{r }
anova(m2,glmFrogs, test="LRT")
```

These p-values are slightly different than the results printed on `summary` which uses a different test (not the "LRT")

```{r }
cat("Summary results:\n")
round(summary(glmFrogs)$coeff,4)
```

We can also use AIC for glm models:

```{r }
AIC(glmFrogs)
```

The AIC is used in the `step` function to do stepwise insertion and deletion of variables, like with regression. 

```{r }
step(glmFrogs, direction="both",trace=0)
```

# Confusion matrices and Classification
The following function `confusion` calculates the entries of the confusion matrix for a given cutoff, given by `thres`. 

```{r }
confusion <- function (y, yhat, thres)
  {
    n <- length(thres)
    conf <- matrix(0,length(thres),ncol=4)
    colnames(conf) <- c("C0","W1","W0","C1")
    rownames(conf)<-round(thres,3)
    for ( i in 1:n)
      {
        C0 <- sum((!y) & (yhat<=thres[i]))
        W1 <- sum((!y) & (yhat>thres[i]))
        W0 <- sum((y) & (yhat<=thres[i])) 
        C1 <- sum((y) & (yhat>thres[i]))
        conf[i,] <- c(C0, W1, W0, C1)
      }
    return(conf)
}
```

We can use this function to calculate the confusion matrix values for a large range of thresholds. The output is a matrix, where each row corresponds to a threshold value, and each columns one of the entries of the 2x2 matrix.

```{r }
thr<-seq(0,1,by=.05)
conf<-confusion(frogs$pres.abs, glmFrogs$fitted.values, thr)
conf
```

Using this output, I can plot the total amount of error, `W1+W0` as a function of the threshold:

```{r }
plot(thr,conf[,"W1"]+conf[,"W0"], xlab="threshold", ylab="W1+W0")
```

# Example of spam data
We will now return to the spam dataset for demonstration. 

```{r }
head(spam7)
```

First we will visualize some of the variables:

```{r }
summary(spam)
par(mfrow=c(3,2))
for (i in 1:5)
  hist(spam[,i],main="",xlab=names(spam)[i],breaks=10000)
par(mfrow=c(1,1))
```

And do a pairs plot:

```{r }
pairs(spam, cex=0.5)
```

The explanatory variables are very skewed, which may be hard for the model. We can consider taking the log-transform of our explanatory variables; however, we will add a small amount before taking the log so we don't run the risk of taking log of 0:

```{r }
s = 0.001  
pairs(~log(crl.tot)+log(dollar+s) +log(bang+s)+log(money+s)+log(n000+s)+log(make+s) + yesno,       data=spam, cex=.5)
```

Based on these plots, we are going to log-transform many of the explanatory variables before putting them into the `glm` model

```{r }
spam.glm <- glm(yesno~ log(crl.tot) + log(dollar+s) + log(bang+s)
                +log(money+s) + log(n000+s) + log(make+s),
                family=binomial, data=spam)
summary(spam.glm)
```

We can compare this to if we had run a standard regression with `lm`:

```{r }
spam.lm <-  lm(as.numeric(yesno=="y")~ log(crl.tot) + log(dollar+s)
               + log(bang+s) +log(money+s) + log(n000+s) + log(make+s)
               ,data=spam)
summary(spam.lm)
```

We can first compare the fitted values from these two models. 

```{r }
par(mfrow=c(1,1))
plot(spam.lm$fitted.values, spam.glm$fitted.values,asp=1)
abline(c(0, 1), col = "red")
```

We can also instead consider the error we make based on the two models. We use our `confusion` function from before to calculate the different type of errors from both models. 

```{r }
v <- seq(0.001,0.999,length=50)
y <- as.numeric(spam$yesno=="y")
glm.conf <- confusion(y,spam.glm$fitted,v)
lm.conf <- confusion(y,spam.lm$fitted,v)
matplot(v,cbind((glm.conf[,"W1"]+glm.conf[,"W0"])/4601, (lm.conf[,"W1"]+lm.conf[,"W0"])/4601), xlab="threshold", ylab="W0+W1",type="b",pch=1)
legend(.8, .4, lty=1:2, col=1:2,c("glm","lm"))
```

We can also compare what would have happened if we hadn't taken the log of our explanatory variables. 

```{r }
spam.glm.nolog <- glm(yesno~ crl.tot + dollar + bang
                +money + n000 + make,
                family=binomial, data=spam)
summary(spam.glm)
summary(spam.glm.nolog)
```

We can see that the AIC is much lower in the one with the log. Similarly with the `lm` models, we have a noticable difference in the R-squared values.

```{r }
spam.lglmFrogs = lm(as.numeric(yesno=="y")~ crl.tot + dollar
                + bang + money + n000 + make, data=spam)
summary(spam.lglmFrogs)
summary(spam.lm)
```

We can use the results of our confusion matrix to calculate a precision/recall curve or a Sensitivity/Specificity plot (also known as a ROC curve):

```{r }
spamGlm.precision<-glm.conf[,"C1"]/(glm.conf[,"C1"]+glm.conf[,"W1"])
spamGlm.recall<-glm.conf[,"C1"]/(glm.conf[,"C1"]+glm.conf[,"W0"])
spamGlm.spec<-glm.conf[,"C0"]/(glm.conf[,"C0"]+glm.conf[,"W1"])
par(mfrow=c(1,2))
plot(x=spamGlm.precision,y=spamGlm.recall,xlab="Precision",ylab="Recall / Sensitivity",type="b",xlim=c(0,1),ylim=c(0,1), main="Precision-Recall")
plot(x=spamGlm.spec,y=spamGlm.recall,ylab="Recall / Sensitivity", xlab="Specificity",type="b",xlim=c(0,1),ylim=c(0,1),main="ROC")
```

Note there is a NaN for precision for the last value of the threshold so there is no point plotted.

We can do the same plots for our linear regression curves:

```{r }
spamlm.precision<-lm.conf[,"C1"]/(lm.conf[,"C1"]+glm.conf[,"W1"])
spamlm.recall<-lm.conf[,"C1"]/(lm.conf[,"C1"]+glm.conf[,"W0"])
spamlm.spec<-lm.conf[,"C0"]/(lm.conf[,"C0"]+glm.conf[,"W1"])

par(mfrow=c(1,2))
matplot(x = cbind(spamGlm.precision,
              spamlm.precision),
        y = cbind(spamGlm.recall,
              spamlm.recall),
        xlab="Precision", ylab="Recall",type="l")
legend(.8, .4, lty=1:2, col=1:2,c("glm","lm"))
matplot(x = cbind(spamGlm.spec,
              spamlm.spec),
        y = cbind(spamGlm.recall,
              spamlm.recall),
        ylab="Recall / Sensitivity", xlab="Specificity",type="l")
legend(.8, .4, lty=1:2, col=1:2,c("glm","lm"))
```

