library(ggplot2)
library(dplyr)
library(tidyverse)
library(knitr)
library(ggthemes)
theme_set(theme_bw())
Steps for Fitting a Model:
Question 1 What is an advantage of RMSE over MSE?
Question 2 When is MSE small? large?
This means that we need to worry about overfitting the training data (doing well on the training data, but poorly on the test data, because our model was too sensitive to outliers)
We have 50 observations on a quantitative response Y and predictor X that we will split 70/30 into training and test data (35/15 observations).
We will fit 3 models with increasing flexibility
This is clearly not linear data, so we would expect a more powerful model to do better.
\(~\)
| model | Train.MSE |
|---|---|
| Linear | 0.677 |
| Quintic | 0.086 |
| Poly | 0.071 |
We build a linear, quintic and 15th degree polynomial model.
| model | Test.MSE |
|---|---|
| Linear | 1.281 |
| Quintic | 0.326 |
| Poly | 1.822 |
Models built on training data are plotted on test data
The 15th degree poly model fits the training data well. But doesn’t do as well on test data.
| model | Train.MSE | Test.MSE |
|---|---|---|
| Linear | 0.677 | 1.281 |
| Quintic | 0.086 | 0.326 |
| Poly | 0.071 | 1.822 |
Question 3 What does the Irreducible error line imply about the distribution of the \(\epsilon\) term in our simulated data?
In addition to the decisions we have to make in Part 1, there are additional tradeoffs that we have to consider when choosing models.
The first tradeoff that the book discusses is the Flexibility-Interpretability tradeoff. I like to think about this as the Accuracy vs Explainable tradeoff. The most powerful models that we have right now are Large Learning Models like GPT. They have billions of parameters (or more), which allows them to model highly complex data. However, they are effectively a black-box model, we don’t really know what’s going on inside them, and require significantly more data (due to the increased parameters). On the other hand, something like OLS regression is easily interpretable, and doesn’t take much data to create an estimated model. However, there are many problems that OLS just wouldn’t work for.
Question 4 The book talks about inference as one reason that you might prefer a less flexible, more interpretable model. Can you think of reasons you might want to use a more interpretable model when making predictions, even though you will be (potentially) less accurate?
Question 5 What types of models tend to have low/high variance? Bias?
Question 6 For the following Research Questions, which choices would you make from the section on “Choices in Model Selection” above? (e.g. would a supervised or unsupervised model work better?). Explain your choices
Which customers are most likely to buy a new product on launch?
How do we identify customer segments (similar customers) based on purchase history?
Which students are likely to pass a course based on grades through week 7?
What is a student’s final (percentage) grade based on grades through week 8?
What will the weather be like tomorrow?
We talked about MSE and RMSE as common error functions, however we could easily come up with others. For example, we could use the Mean Error (no squaring, this is a bad idea, but a good thought experiment), the Mean Absolute Error (MAE), the Mean Quartic Error (4th power), or the Max error.
Question 7 Assuming you are fitting a line, how do you think changing the error functions would effect your line with respect to accuracy on the training set and accuracy on outliers in the training set?
Question 8 Using the code below as a starting place, explore how the choice of the error function effects the model. Do your predictions match your expectations? Why or Why not? Were you surprised by any of this? In addition, what is \(f(x)\) and what is \(\epsilon\) (the irreducible error)?
Note: You should not need to change the following sections:
You will need to change the defining the objective section (uncomment the line you want to use).
You may want to change the plotting section (e.g. to plot multiple outputs)
Bonus Question 1 If time permits, try generating your own error (objective) functions.
library(CVXR)
# Generate Data
set.seed(30)
X <- runif(50, 0,10)
Y <- 2*X+ rnorm(50, 0, 2)
my_df <- data.frame(X,Y)
# defining variables to be tuned during optimisation
coefficient <- Variable(1)
intercept <- Variable(1)
# defining the objective i.e. minimizing the sum af absolute differences (MAE)
# This is the main part that you should change
#objective <- Minimize(sum(abs(my_df$Y - (my_df$X * coefficient) - intercept))) #MAE
objective <- Minimize(sum((my_df$Y - (my_df$X * coefficient) - intercept)**2)) #MSE
#objective <- Minimize(sum((my_df$Y - (my_df$X * coefficient) - intercept)**4)) #MQE
#objective <- Minimize(max(abs(my_df$Y - (my_df$X * coefficient) - intercept))) #Max Error
# optimization
problem <- Problem(objective)
result <- solve(problem)
# result
result$status #This should say "optimal"
## [1] "optimal"
result_coefficient <- result$getValue(coefficient)
result_intercept <- result$getValue(intercept)
# Plot the data
Yhat <-(my_df$X * result_coefficient) + result_intercept
my_df2 <- data.frame(X,Y,Yhat)
ggplot(my_df2, aes(x=X))+
geom_point(aes(y=Y),color="black")+
geom_line(aes(y=Yhat),color="blue")
If time permits, I would like you to explore this section. The following are useful papers if you would like to read about them in the future. The summary is that different fairness definitions (that are all reasonable) are mutually incompatible in most cases.
This section will focus on the paper “On Fairness and Calibration” which shows “that calibration is compatible only with a single error constraint (i.e. equal false-negatives rates across groups), and show that any algorithm that satisfies this relaxation is no better than randomizing a percentage of predictions for an existing classifier”
This is a simpler “proof” of the Impossibility of Equalized Odds with Calibration from “Inherent Trade-Offs in the Fair Determination of Risk Scores”. In other words, it is impossible to balance the standard rates when the underlying groups are different (except in trivial cases). For example, Looking at (b) and (c) we can see that h1 and h2 cannot be balanced under both False Positive and False Negative rates
The following is taken directly from “Inherent Trade-Offs in the Fair Determination of Risk Scores”. (with some formatting changes to make it more legible in Rstudio)
“To take one simple example, suppose we want to determine the risk that a person is a carrier for a disease X, and suppose that a higher fraction of women than men are carriers. Then our results imply that in any test designed to estimate the probability that someone is a carrier of X, at least one of the following undesirable properties must hold:
The point is that this trade-off among (a), (b), and (c) is not a fact about medicine; it is simply a fact about risk estimates when the base rates differ between two groups.”
In other words, if the underlying groups aren’t the same, we can’t balance the rates (True/False Positive/Negative) between the groups (except in the trivial cases). Below is a (made up example).
# Load and Clean the base data
library(mlbench)
## Warning: package 'mlbench' was built under R version 4.4.2
data(BreastCancer)
BC2<-BreastCancer[,-1] #Drop Id
BC2<-drop_na(BC2) #Drop NA
BC2<-lapply(BC2, function(x) type.convert(as.character(x), as.is = TRUE)) #Drop Levels
BC2$Class<-BC2$Class=="benign" #Create Boolean Outcome
# Make it fake by adding a "ProtClass" Variable to demonstrate that we can't make all 4 rates equal
set.seed(30)
ProtClass=1
for (i in 2:length(BC2$Class)){
temp=BC2$Class[i]
temp_prob=runif(1)
if (temp==1){
if(temp_prob>.45) {#more likely to have cancer given everything else
ProtClass=c(ProtClass,1)
}
else{
ProtClass=c(ProtClass,0)
}
}
else{
if(temp_prob>.5){#same likelihood to not
ProtClass=c(ProtClass,1)
}
else{
ProtClass=c(ProtClass,0)
}
}
}
BC2$ProtClass<-ProtClass
# Fit the Log Model
pred_model<-glm(Class ~ .,data=BC2,family=binomial)
pred_out<-predict(pred_model,type="response")
Y_Yhat<-data.frame(BC2$Class,pred_out,BC2$ProtClass)
ggplot(Y_Yhat, aes(x=BC2.Class,y=pred_out,color=BC2.ProtClass))+
geom_point()
We can see that this data is not easily separable by T/F or by
Class.
\(~\) Below you can see how I generated the rates and accuracy for all whole percentage thresholds. You are free to play with anything below if there are other values you’d like to see.
# Calculate all Rates Based on 1% Thresholds (for Positivity) i.e. all predicted values greater than threshold are set to 1)
x=1:100/100
Rates_singular<-function(threshold,data){
temp=data$pred_out>=threshold
TP=sum(data$BC2.Class==1&temp==1)
FP=sum(data$BC2.Class==0&temp==1)
FN=sum(data$BC2.Class==1&temp==0)
TN=sum(data$BC2.Class==0&temp==0)
Acc=(TP+TN)/(TP+FP+FN+TN)
return(c(TP/(TP+FN),FP/(TN+FP),FN/(TP+FN),TN/(TN+FP),Acc))
}
Rates<-function(threshold, data){
temp0=data[data$BC2.ProtClass==0,]
temp1=data[data$BC2.ProtClass==1,]
return(c(Rates_singular(threshold,temp0),Rates_singular(threshold,temp1)))
}
StatRates<-data.frame(x=0,TP0=0,FP0=0,FN0=0,TN0=0,Acc0=0,TP1=0,FP1=0,FN1=0,TN1=0,Acc1=0)
StatRates[1,]<-c(0,Rates(0,Y_Yhat))
for(i in 1:length(x)){
StatRates=rbind(StatRates,c(x[i],Rates(x[i],Y_Yhat)))
}
Plot the data. This is the part that you should play with. Other than the trivial cases where we aren’t very accurate (e.g. assume everyone does/doesn’t have cancer), How would you balance the rates to be fair? You should look at the various rates you could plot, as well as the limits of the graph.
ggplot(StatRates, aes(x=x))+
geom_line(aes(y=TP0),color="blue")+
geom_line(aes(y=TP1),color="red")+
geom_line(aes(y=Acc0),color="black",size=1)+
geom_line(aes(y=Acc1),color="green",size=1)+
xlim(.4,.6)+
ylim(.9,1)
## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
Finding the “best” value by Rate we wish to optimize. E.g. If we wanted to maximize accuracy for Group 0, we’d want to use a 58% threshold, while Group 1 would want an 83% threshold.
for (col in colnames(StatRates)){
print(paste(col, which.min(StatRates[,col]), which.max(StatRates[,col])))
}
## [1] "x 1 101"
## [1] "TP0 101 1"
## [1] "FP0 84 1"
## [1] "FN0 1 101"
## [1] "TN0 1 84"
## [1] "Acc0 101 58"
## [1] "TP1 101 1"
## [1] "FP1 98 1"
## [1] "FN1 1 101"
## [1] "TN1 1 98"
## [1] "Acc1 101 83"
In today’s lecture and lab we covered the following: