I am trying to estimate a logistic regression, using the 10-fold cross-validation.
#import libraries
library(car); library(caret); library(e1071); library(verification)
#data import and preparation
data(Chile)
chile <- na.omit(Chile) #remove "na's"
chile <- chile[chile$vote == "Y" | chile$vote == "N" , ] #only "Y" and "N" required
chile$vote <- factor(chile$vote) #required to remove unwanted levels
chile$income <- factor(chile$income) # treat income as a factor
Goal is to estimate a glm - model that predicts to outcome of vote "Y" or "N" depended on relevant explanatory variables and, based on the final model, compute a confusion matrix and ROC curve to grasp the models behaviour for different threshold levels.
Model selection leads to:
res.chileIII <- glm(vote ~
sex +
education +
statusquo ,
family = binomial(),
data = chile)
#prediction
chile.pred <- predict.glm(res.chileIII, type = "response")
generates:
> head(chile.pred)
1 2 3 4 5 6
0.974317861 0.008376988 0.992720134 0.095014139 0.040348115 0.090947144
to compare the observed with estimation:
chile.v <- ifelse(chile$vote == "Y", 1, 0) #to compare the two arrays
chile.predt <- function(t) ifelse(chile.pred > t , 1,0) #t is the threshold for which the confusion matrix shall be computed
confusion matrix for t = 0.3:
confusionMatrix(chile.predt(0.3), chile.v)
> confusionMatrix(chile.predt(0.3), chile.v)
Confusion Matrix and Statistics
Reference
Prediction 0 1
0 773 44
1 94 792
Accuracy : 0.919
95% CI : (0.905, 0.9315)
No Information Rate : 0.5091
P-Value [Acc > NIR] : < 2.2e-16
and the Roc-curve:
roc.plot(chile.v, chile.pred)
which seems as a reasonable model.
Now instead of using the "normal" predict.glm() function I want to test out the performance difference to a 10-fold cross-validation estimation.
tc <- trainControl("cv", 10, savePredictions=T) #"cv" = cross-validation, 10-fold
fit <- train(chile$vote ~ chile$sex +
chile$education +
chile$statusquo ,
data = chile ,
method = "glm" ,
family = binomial ,
trControl = tc)
> summary(fit)$coef
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.0152702 0.1889646 5.372805 7.752101e-08
`chile$sexM` -0.5742442 0.2022308 -2.839549 4.517738e-03
`chile$educationPS` -1.1074079 0.2914253 -3.799971 1.447128e-04
`chile$educationS` -0.6827546 0.2217459 -3.078996 2.076993e-03
`chile$statusquo` 3.1689305 0.1447911 21.886224 3.514468e-106
all parameters significant.
fitpred <- ifelse(fit$pred$pred == "Y", 1, 0) #to compare with chile.v
> confusionMatrix(fitpred,chile.v)
Confusion Matrix and Statistics
Reference
Prediction 0 1
0 445 429
1 422 407
Accuracy : 0.5003
95% CI : (0.4763, 0.5243)
No Information Rate : 0.5091
P-Value [Acc > NIR] : 0.7738
which is obviously very different from the previous confusion matrix. My expectation was that the cross validated results should not perform much worse then the first model. However the results show something else.
My assumption is that there is a mistake with the settings of the train() parameters but I can't figure it out what it is.
I would really appreciate some help, thank you in advance.
You are trying to get an idea of the in-sample fit using a confusion matrix. Your first approach using the glm() function is fine.
The problem with the second approach using train() lies in the returned object. You are trying to extract the in-sample fitted values from it by fit$pred$pred. However, fit$pred does not contain the fitted values that are aligned to chile.v or chile$vote. It contains the observations and fitted values of the different (10) folds:
head(fit$pred)
# pred obs rowIndex parameter Resample
# 1 N N 2 none Fold01
# 2 Y Y 20 none Fold01
# 3 Y Y 28 none Fold01
# 4 N N 38 none Fold01
# 5 N N 55 none Fold01
# 6 N N 66 none Fold01
tail(fit$pred)
# pred obs rowIndex parameter Resample
# 1698 Y Y 1592 none Fold10
# 1699 Y N 1594 none Fold10
# 1700 N N 1621 none Fold10
# 1701 N N 1656 none Fold10
# 1702 N N 1671 none Fold10
# 1703 Y Y 1689 none Fold10
Due to the randomness of the folds, and because you are predicting 0 or 1, you get an accuracy of roughly 50%.
The in-sample fitted values you are looking for are in fit$finalModel$fitted.values. Using those:
fitpred <- fit$finalModel$fitted.values
fitpredt <- function(t) ifelse(fitpred > t , 1,0)
confusionMatrix(fitpredt(0.3),chile.v)
# Confusion Matrix and Statistics
#
# Reference
# Prediction 0 1
# 0 773 44
# 1 94 792
#
# Accuracy : 0.919
# 95% CI : (0.905, 0.9315)
# No Information Rate : 0.5091
# P-Value [Acc > NIR] : < 2.2e-16
#
# Kappa : 0.8381
# Mcnemar's Test P-Value : 3.031e-05
#
# Sensitivity : 0.8916
# Specificity : 0.9474
# Pos Pred Value : 0.9461
# Neg Pred Value : 0.8939
# Prevalence : 0.5091
# Detection Rate : 0.4539
# Detection Prevalence : 0.4797
# Balanced Accuracy : 0.9195
#
# 'Positive' Class : 0
Now the accuracy is around the expected value. Setting the threshold to 0.5 yields about the same accuracy as the estimate from the 10-fold cross-validation:
confusionMatrix(fitpredt(0.5),chile.v)
# Confusion Matrix and Statistics
#
# Reference
# Prediction 0 1
# 0 809 64
# 1 58 772
#
# Accuracy : 0.9284
# 95% CI : (0.9151, 0.9402)
# [rest of the output omitted]
fit
# Generalized Linear Model
#
# 1703 samples
# 7 predictors
# 2 classes: 'N', 'Y'
#
# No pre-processing
# Resampling: Cross-Validated (10 fold)
#
# Summary of sample sizes: 1533, 1532, 1532, 1533, 1532, 1533, ...
#
# Resampling results
#
# Accuracy Kappa Accuracy SD Kappa SD
# 0.927 0.854 0.0134 0.0267
Additionally, regarding your expectation "that the cross validated results should not perform much worse than the first model", please check summary(res.chileIII) and summary(fit). The fitted models and coefficients are exactly the same so they will give the same results.