K-fold cross-validation is a re-sampling technique that we use to assess the performance of a machine learning model. It works by separating the data into k groups of equal size. Each fold serves as a validation set once, while the remaining k folds serve as the training set. Then, for each fold, the model is trained on the remaining folds, and then tested on that separated fold while we store the performance metrics like test error, RMSE, etc. After k iterations, we take the average of the performance metrics to get a more reliable estimate of model performance.
K-fold cross-validation is more reliable than the validation set approach because it uses multiple train-test splits, which reduces the variance and provides a better estimate of model performance. Validation set approach only uses 1 split, which could lead to biased results, especially if the chosen validation set is not representative of the entire dataset.
On the other hand, K-fold cross-validation is much more computationally expensive than the validation set approach. The higher number of k, the more computational power is required and the longer it would take to implement. Also, k-fold cross-validation ensures that all data points are used for both training and testing, which makes it a better choice for small datasets where data efficiency is crucial. The validation set approach is faster and easier to implement, which makes it useful if computational resources are limited.
LOOCV is basically the k-fold cross-validation where our k = our number of observations so that each fold is 1 observation, and it uses all of the data except 1 point which is for validation. This leads to very minimal bias, but a very high computational cost. While LOOCV provides a more exhaustive use of the data, it can result in high variance in predictions and it is sort of impractical for large datasets.
Load ISLR2 and Default dataset
library(ISLR2)
library(caret)## Loading required package: ggplot2## Loading required package: latticedata('Default')set.seed(42)
log_model_5a <- glm(default ~ income + balance,
data = Default, family = 'binomial')
summary(log_model_5a)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.154e+01 4.348e-01 -26.545 < 2e-16 ***
## income 2.081e-05 4.985e-06 4.174 2.99e-05 ***
## balance 5.647e-03 2.274e-04 24.836 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1579.0 on 9997 degrees of freedom
## AIC: 1585
##
## Number of Fisher Scoring iterations: 8train_index_default <- createDataPartition(Default$default, p = 0.7, list = FALSE)
train_data_default <- Default[train_index_default, ]
test_data_default <- Default[-train_index_default, ]log_model_5a_ii <- glm(default ~ income + balance,
data = train_data_default, family = 'binomial')pred_probs_5a_ii <- predict(log_model_5a_ii, newdata = test_data_default, type = 'response')
predictions_5a_ii <- ifelse(pred_probs_5a_ii > 0.5, 'Yes', 'No')mean(predictions_5a_ii != test_data_default$default)## [1] 0.02734245The test error rate for the validation set approach is equal to 2.73% using seed 42.
set.seed(123)
#New split with different seed
train_index_default2 <- createDataPartition(Default$default, p = 0.7, list = FALSE)
train_data_default2 <- Default[train_index_default2, ]
test_data_default2 <- Default[-train_index_default2, ]
#New model
log_model_5a_ii_2 <- glm(default ~ income + balance,
data = train_data_default2, family = 'binomial')
#New predictions
pred_probs_5a_ii_2 <- predict(log_model_5a_ii_2, newdata = test_data_default2,
type = 'response')
predictions_5a_ii_2 <- ifelse(pred_probs_5a_ii_2 > 0.5, 'Yes', 'No')
mean(predictions_5a_ii_2 != test_data_default2$default)## [1] 0.02534178The test error rate for the validation set approach is equal to 2.53% using seed 123.
set.seed(1)
#New split with different seed
train_index_default3 <- createDataPartition(Default$default, p = 0.7, list = FALSE)
train_data_default3 <- Default[train_index_default3, ]
test_data_default3 <- Default[-train_index_default3, ]
#New model
log_model_5a_ii_3 <- glm(default ~ income + balance,
data = train_data_default3, family = 'binomial')
#New predictions
pred_probs_5a_ii_3 <- predict(log_model_5a_ii_3, newdata = test_data_default3,
type = 'response')
predictions_5a_ii_3 <- ifelse(pred_probs_5a_ii_3 > 0.5, 'Yes', 'No')
mean(predictions_5a_ii_3 != test_data_default3$default)## [1] 0.02600867The test error rate for the validation set approach is equal to 2.60% using seed 1.
Using 3 different splits showed 3 different results as expected. The main issue is that in our initial split we got a test error of 2.73%, and on our following 2 splits we got 2.53% and 2.60%. This shows a small variability in the results but it is safe to say that our model is not overly sensitive to particular train-validation split and our error rate is likely a reliable estimate.
set.seed(42)
#New split with different seed
train_index_default_5d <- createDataPartition(Default$default, p = 0.7, list = FALSE)
train_data_default_5d <- Default[train_index_default_5d, ]
test_data_default_5d <- Default[-train_index_default_5d, ]
log_model_5d <- glm(default ~ income + balance + student,
data = train_data_default_5d, family = 'binomial')
summary(log_model_5d)##
## Call:
## glm(formula = default ~ income + balance + student, family = "binomial",
## data = train_data_default_5d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.113e+01 5.923e-01 -18.786 <2e-16 ***
## income 5.857e-06 9.767e-06 0.600 0.5487
## balance 5.837e-03 2.831e-04 20.618 <2e-16 ***
## studentYes -5.633e-01 2.854e-01 -1.974 0.0484 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2050.6 on 7000 degrees of freedom
## Residual deviance: 1076.9 on 6997 degrees of freedom
## AIC: 1084.9
##
## Number of Fisher Scoring iterations: 8#New predictions
pred_probs_5d <- predict(log_model_5d, newdata = test_data_default_5d,
type = 'response')
predictions_5d <- ifelse(pred_probs_5d > 0.5, 'Yes', 'No')
mean(predictions_5d != test_data_default_5d$default)## [1] 0.02834278The test error rate for the validation set approach is equal to 2.83% with seed 42.
set.seed(123)
#New split with different seed
train_index_default_5d_2 <- createDataPartition(Default$default, p = 0.7, list = FALSE)
train_data_default_5d_2 <- Default[train_index_default_5d_2, ]
test_data_default5d_2 <- Default[-train_index_default_5d_2, ]
log_model_5d_2 <- glm(default ~ income + balance + student,
data = train_data_default_5d, family = 'binomial')
#New predictions
pred_probs_5d_2 <- predict(log_model_5d_2, newdata = test_data_default5d_2,
type = 'response')
predictions_5d_2 <- ifelse(pred_probs_5d_2 > 0.5, 'Yes', 'No')
mean(predictions_5d_2 != test_data_default5d_2$default)## [1] 0.02567523The test error rate for the validation set approach is equal to 2.57% with seed 123.
After testing the model with the added student variable, we can see that in both splits the results were very similar, arguably a bit worse. We know that when adding ‘student’ our ‘income’ variable becomes insignificant so in terms of performance, we didn’t expect it to be better since we added a variable that made another one have basically no effect on ‘default’.
In summary, adding ‘student’ does not lead to a reduction in test error.
set.seed(42)
log_model_6 <- glm(default ~ income + balance,
data = Default, family = 'binomial')
summary(log_model_6)##
## Call:
## glm(formula = default ~ income + balance, family = "binomial",
## data = Default)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.154e+01 4.348e-01 -26.545 < 2e-16 ***
## income 2.081e-05 4.985e-06 4.174 2.99e-05 ***
## balance 5.647e-03 2.274e-04 24.836 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2920.6 on 9999 degrees of freedom
## Residual deviance: 1579.0 on 9997 degrees of freedom
## AIC: 1585
##
## Number of Fisher Scoring iterations: 8Estimated standard error for income:
0.00002081
Estimated standard error for balance:
0.005647
This means that for an increase of 1 unit in income, the log-odds of default being ‘yes’ go increase by 0.00002081
This means that for an increase of 1 unit in balance, the log-odds of default being ‘yes’ increase by 0.005647
library(boot)##
## Attaching package: 'boot'## The following object is masked from 'package:lattice':
##
## melanoma#define our boot.fn function
boot.fn <- function(data, index) {
model <- glm(default ~ income + balance,
data = data[index, ],
family = 'binomial')
#return the coefficients for income and balance
return(coef(model)[c('income', 'balance')])
} set.seed(42)
boot(Default, boot.fn, R = 1000)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Default, statistic = boot.fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 2.080898e-05 2.737444e-08 5.073444e-06
## t2* 5.647103e-03 1.176249e-05 2.299133e-04Standard errors from glm() Function:
income = 4.985e-06
balance = 2.274e-04
Standard errors from Bootstrap:
income = 5.073444e-06
balance = 2.299133e-04
Standard errors for both income
and balance are very close to each
other in the glm() function and Bootstrap.
This indicates consistency and reliability between the two approaches. In this case, the bootstrap method does not provide significantly different results, which means the logistic regression model is well specified since bootstrap is more robust and shows extremely similar results.
data(Boston)summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio lstat
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 1.73
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.: 6.95
## Median : 5.000 Median :330.0 Median :19.05 Median :11.36
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :12.65
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:16.95
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :37.97
## medv
## Min. : 5.00
## 1st Qu.:17.02
## Median :21.20
## Mean :22.53
## 3rd Qu.:25.00
## Max. :50.00mu <- mean(Boston$medv)
print(mu)## [1] 22.53281Hint: We can compute the standard error of the sample mean by dividing the sample standard deviation by the square root of the number of observations.
std_error_mu <- sd(Boston$medv) / sqrt(length(Boston$medv))
print(std_error_mu)## [1] 0.4088611boot.b <- function(data, index) {
return(mean(data[index]))
}
set.seed(42)
boot_results <- boot(Boston$medv, boot.b, R = 1000)
boot_se <- sd(boot_results$t) #extract bootstrap standard error
boot_results #print bootstrap results##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.b, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 22.53281 0.02671186 0.4009216The bootstrap standard error does not differ almost at all from the formula based standard error.
Bootstrap std.error = 0.4009216
Formula std.error = 0.4088611
t.test(Boston$medv)##
## One Sample t-test
##
## data: Boston$medv
## t = 55.111, df = 505, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 21.72953 23.33608
## sample estimates:
## mean of x
## 22.53281confidence_interval <- c(mu - 2 * boot_se, mu + 2 * boot_se)
confidence_interval## [1] 21.73096 23.33465Just like our standard errors, confidence intervals are very close in both our bootstrap estimate and the one from t.test(Boston$medv).
t.test confidence intervals: 21.72952, 23.33608
bootstrap confidence intervals: 21.73096, 23.33456
mu_med <- median(Boston$medv)
mu_med## [1] 21.2boot.m <- function(data, index) {
return(median(data[index]))
}
set.seed(42)
boot_results_median <- boot(Boston$medv, boot.m, R = 1000)
boot_se_median <- sd(boot_results_median$t) #extract bootstrap standard error
boot_results_median #print results of bootstrap##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.m, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 21.2 0.0106 0.3661785We get a standard error of 0.361785 for ^µmed.
This indicates how much the median estimate varies across different resamples
mu_01 <- quantile(Boston$medv, 0.1)
mu_01## 10%
## 12.75boot.fn.percentile <- function(data, index) {
return(quantile(data[index], 0.1))
}
set.seed(42)
boot_results_01 <- boot(Boston$medv, boot.fn.percentile, R = 1000)
boot_se_01 <- sd(boot_results_01$t)
boot_results_01##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = Boston$medv, statistic = boot.fn.percentile, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 12.75 0.0215 0.4948966The standard error for ˆ µ0.1 is 0.4948966.
This indicates that 10th percentile estimate varies that much across different resamples.