Checking Shape

When we build logistic regression models, we do make a shape assumption that the relationship between \(X\) and the log odds of the penguin being male is a line. To check this, we create a special scatter plot called an empirical log odds plot. To teach R the function it needs to make this plot, create an R chunk, copy the function below, paste it into an R chunk in your Markdown, and press play.

logOddsPlot <- function(x, y, xname, formulahere){

  if(class(y)=="factor"){
    baseline = levels(y)[1]
    y <- ifelse(y == baseline, 0, 1)
  }
  
  sort = order(x)
  x = x[sort]
  y = y[sort]
  a = seq(1, length(x), by=.05*length(y))
  b = c(a[-1] - 1, length(x))
  
  prob = xmean = ns = rep(0, length(a)) # ns is for CIs
  for (i in 1:length(a)){
    range = (a[i]):(b[i])
    prob[i] = mean(y[range])
    xmean[i] = mean(x[range])
  }
  
  extreme = (prob == 1 | prob == 0)
  prob[prob == 0] = .01
  prob[prob == 1] = .99
  
  g = log(prob/(1-prob))
  
  dataHere <- data.frame("x" = b, "LogOdds" = g)
  
  suppressMessages(library(ggplot2))
  
  ggplot(dataHere, aes(x =xmean, y = LogOdds)) + geom_point() + geom_smooth(formula = formulahere, method = "lm", se = FALSE) + labs(x = xname, y = "Log Odds")
}

To use the function, you will need to specify the inputs, meaning the pieces of information the computer needs in order to run the code.

• x: the column containing the x variable; example: penguins[, "flipper_length_mm"]
• y: the column containing the y variable; example: penguins[, "sex"]
• xname: the label you want on the x axis, example: xname = "Bill Depth in mm". - formula: This is the shape of the relationship you want to plot. If you want a line, this is y ~ x.

Here is some of the code to create the plot. You will need to fill in the gaps.

logOddsPlot(x= , y= , xname = "Bill Depth (mm)", formulahere = y ~ x)

Typically, we would need to check the shape condition for every numeric feature. For today, you may proceed as though the shape conditions is met. When you work on your final project, you will need to check shape if you choose to use logistic regression or linear regression. Be aware that you may find when you check shape that you need to transform X in order to make sure the shape of the relationship is accurately captured in the model.

Logistic Regression in R

The example logistic regression model only uses one feature, but logistic regression models can use a combination of categorical and numeric features, just like every other model in our course.

To build a logistic regression model in R using all of the features in our data set, we use the following code:

m1 <- glm( sex ~ . , data = penguins, family = "binomial")

If we want to see the coefficients from the model, we can use the same function we would use for linear regression:

summary(m1)

These coefficients help us use the model for both association and prediction tasks. They describe the relationship between the features and the log odds that a penguin is male.

Our goal for today, though, is prediction. This means we are not interested in interpreting the coefficients. Instead, we want to use the model to predict whether or not a penguin is male based on the features.

The first step to doing this is to use our trained model to find \(\hat{\pi}_i\) for each penguin. In order to find the predicted probability that each penguin in the data set is male, we use the predict function (again, standard in our course):

probabilities <- predict(m1, type = "response")

The only new part of our code is type = "response". This is important because in logistic regression we might want the predicted log odds or the predicted probabilities to be produced with this function. type = response tells R we can the predicted probabilities.

Now, our goal for today is prediction, which means we need to convert those probabilities \(\hat{\pi}_i\) into predictions \(\hat{y}_i\) for each row \(i\) in the data set. This means we need to decide that probability \(\hat{\pi}_i = P(\hat{y}_i = 1)\) is large enough that we want to assign \(\hat{y}_i = 1\).

This probability that is “large enough” is called a threshold. This is a value such that if the probability that a penguin \(i\) is male (\(\hat{\pi}_i\)) is above the threshold, we predict \(\hat{y}_i = male\). If \(\hat{\pi}_i\) is less than or equal to the threshold, we predict \(\hat{y}_i = female\).

Once we know our threshold, we can use R to convert all the probabilities in the vector probabilities into predictions.

YHat_logistic <- ifelse( probabilities > .5, "male", "female")

At this point, we have vector YHat_logistic that holds a predicted value of \(Y\) for all rows in the data set. The next step is to determine how accurately these predictions reflect the true values in the data set!

Assessing Accuracy

There are a lot of tools we use to assess the predictive accuracy of a classification task, including the accuracy, the classification error rate (CER), the true positive rate (TPR), and the true negative rate (TNR). It is common to consider TPR and TNR when we think about classification tasks, so we will focus on those for today.

TPR is the true positive rate. This phrase makes sense when our response variable is “yes” and “no”, but it’s trickier when we have data like ours when \(Y\) is “male” or “female”. “Positive” always refers to whatever level of \(Y\) is classified as a 1. This means that for our data, the TPR is the percent of the male penguins in the original data set that we correctly predict as male using our logistic regression model with a specific threshold.

\[TPR = \frac{ \sum_{i=1}^n (y_i = male ~ and ~ \hat{y}_i = male)}{\sum_{i=1}^n (y_i = male)}\]

To compute the TPR, we can use the following function:

compute_TPR <- function( truth, predictions, isone){
  
  # Determine which rows in the true data have y = 1
  true1 <- which( truth == isone)
  
  # Determine which rows in the predictions have hat y = 1
  pred1 <- which( predictions == isone)
  
  # Find how many agree
  numerator <- intersect( true1, pred1)
  
  # Compute the TPR
  TPR <- length(numerator)/length(true1)
  
  TPR
  
  
}

This function requires 3 inputs:

• truth: What is the actual \(Y\) variable we trying to predict?
• predictions: what are the predicted values \(\hat{Y}\) from our model?
• isone: what value of the response variable are we computing for? In this case, we want the percent of the male penguins that are correctly estimated, so isone = "male".

We use the function like so:

compute_TPR( truth = penguins[, "sex"], predictions = YHat_logistic, isone = "male")

In Question 5, we will notice something very odd. Our TPR is very high for very small choice of threshold. Why is that??

Well, the smaller the threshold, the easier we make it to declare that \(\hat{y}_i= male\). This means that the probability that a penguin that actually is male is predicted to be male is very high, and this is what the TPR measures!

However, this means that it is also very likely that we assign the value \(\hat{y}_i=male\) when the true value of \(y_i\) is female. This is why in addition to the TPR, we discuss the TNR. In our case, the TNR is the percent of the female penguins in the original data set that we correctly predict as female using our logistic regression model with a specific threshold.

\[TNR = \frac{ \sum_{i=1}^n (y_i = female ~ and ~ \hat{y}_i = female)}{\sum_{i=1}^n (y_i = female)}\]

As we can see, though our TPR is high with a threshold of .001, our TNR is not. We typically need a balance between the TPR and the TNR - we want to predict well for both male and female penguins, so we want both values to be relatively high!

Because of this, there are a few different tools we can use that assess how well our model is doing at prediction that take both TPR and TNR into account. One such tool is the geometric mean of TPR and TNR:

\[\sqrt{TPR \times TNR}\]

The quantity above tends to be closer to 1 when both TPR and TNR are high, meaning both \(Y_i =0\) and \(Y_i = 1\) are being predicted well. The quantity is closer to 0 when at least one of \(Y_i =0\) and/or \(Y_i = 1\) is NOT being predicted well. This means the model is doing better at prediction the closer the geometric mean of TPR and TNR gets to 1.

Another score that is sometimes used in classification is the F1 Score, which is defined as:

\[F1 = \frac{2(TPR)}{2(TPR) + FPR + FNR}\]

where

• \(TPR\) = the proportion of male penguins that were correct predicted as males
• \(FPR\) = the proportion of female penguins that were incorrectly predicted as males.

The highest possible value for the F1 score is 1, while the lowest is 0, just like the geometric mean of TPR and TNR. This means models with higher F1 score have higher predictive accuracy.

This means we need TPR, FNR, FPR…and we really need a confusion matrix.

knitr::kable( table( YHat_logistic, penguins$sex) ,col.names = c("True Female", "True Male"))

Note: In this output, the TPR = sensitivity, the TNR = specificity.

One thing to note is that so far, all of the values we have computed tell about the in-sample predictive ability of our model. If we want to estimate how well our model will predict for test data, we need to use a cross validation technique like LOOCV or 10 fold-CV, or we need test data.

My folks using classification on your project, our goal is prediction, so you will need to use CV to estimate your predictive accuracy!!! This means your code will be a little different from what we do in this lab!!

Predicting with a tree

To obtain the predicted values from the tree, we use the following:

YHat_tree <- predict(m3, type = "class")