R - cv.glm$delta vs. semi-manual calculation of mean squared error I went to https://www.lendingclub.com/info/download-data.action?file=LoanStats.csv and downloaded the file LoanStats_2017Q4.csv. I sampled 2% of the data and ran it through some cleaning (177 lines of code that removed qualitative variables, set the string binaries to a numeric 1/0, turned strings into numbers, removed obvious linear combinations, and omitted records with "na"s in them) to go from a 118648x145 data file to 1818x71. I'm attempting to predict "int_rate" based on the other variables. Per request here is the cleaning code I used to generate my 1818x71 dataset from the full size dataset. setwd("[Your_storage_location") project_data <- read.csv("LoanStats_2017Q4_APM_Project.csv", header = TRUE, stringsAsFactors = FALSE) #bring in the data dim(project_data) #118648x145 set.seed(42) reduce_to_make_processing_easier <- sample(x = c(TRUE,FALSE), size = (dim(project_data)[1]), rep=TRUE, prob = c(.02,.98)) project_data_reduced1 <- project_data[reduce_to_make_processing_easier,] dim(project_data_reduced1) #2441x145 del_cols_1 <- c("id","member_id","grade","sub_grade","home_ownership", "policy_code","collection_recovery_fee","purpose", "loan_status","verification_status","home_ownership", "emp_title","issue_d","url","desc","title","zip_code", "addr_state","earliest_cr_line","last_pymnt_d", "next_pymnt_d","last_credit_pull_d","hardship_start_date", "hardship_end_date","payment_plan_start_date", "debt_settlement_flag_date","settlement_date", "mths_since_last_delinq","mths_since_last_record", "annual_inc_joint", "dti_joint", "verification_status_joint", "revol_bal_joint", "sec_app_earliest_cr_line", "sec_app_inq_last_6mths", "sec_app_mort_acc", "sec_app_open_acc", "sec_app_revol_util", "sec_app_open_act_il", "sec_app_num_rev_accts", "sec_app_chargeoff_within_12_mths", "sec_app_collections_12_mths_ex_med", "sec_app_mths_since_last_major_derog", "hardship_type", "hardship_reason", "hardship_status", "deferral_term", "hardship_amount", "hardship_start_date", "hardship_end_date", "payment_plan_start_date", "hardship_length", "hardship_dpd", "hardship_loan_status", "orig_projected_additional_accrued_interest", "hardship_payoff_balance_amount", "hardship_last_payment_amount", "debt_settlement_flag_date", "settlement_status", "settlement_date", "settlement_amount", "settlement_percentage", "settlement_term", "mths_since_last_major_derog", "bc_open_util", "bc_util", "mths_since_recent_revol_delinq", "mths_since_recent_bc_dlq", "emp_length") project_data_reduced <- project_data_reduced1[,!(names(project_data_reduced1) %in% del_cols_1)] dim(project_data_reduced) #2441x83 #reset binary columns to a 1/0 dsf <- which(colnames(project_data_reduced)=="debt_settlement_flag") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,dsf] == "N") { project_data_reduced[row,dsf] = 0 } else { project_data_reduced[row,dsf] = 1 } } project_data_reduced[,dsf] <- sapply(project_data_reduced[,dsf], as.integer) dm <- which(colnames(project_data_reduced)=="disbursement_method") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,dm] == "DirectPay") { project_data_reduced[row,dm] = 0 } else { project_data_reduced[row,dm] = 1 } } project_data_reduced[,dm] <- sapply(project_data_reduced[,dm], as.integer) hf <- which(colnames(project_data_reduced)=="hardship_flag") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,hf] == "N") { project_data_reduced[row,hf] = 0 } else { project_data_reduced[row,hf] = 1 } } project_data_reduced[,hf] <- sapply(project_data_reduced[,hf], as.integer) apt <- which(colnames(project_data_reduced)=="application_type") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,apt] == "Joint app") { project_data_reduced[row,apt] = 0 } else { project_data_reduced[row,apt] = 1 } } project_data_reduced[,apt] <- sapply(project_data_reduced[,apt], as.integer) ils <- which(colnames(project_data_reduced)=="initial_list_status") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,ils] == "w") { project_data_reduced[row,ils] = 0 } else { project_data_reduced[row,ils] = 1 } } project_data_reduced[,ils] <- sapply(project_data_reduced[,ils], as.integer) pp <- which(colnames(project_data_reduced)=="pymnt_plan") for (row in 1:dim(project_data_reduced)[1]) { if (project_data_reduced[row,pp] == "n") { project_data_reduced[row,pp] = 0 } else { project_data_reduced[row,pp] = 1 } } project_data_reduced[,pp] <- sapply(project_data_reduced[,pp], as.integer) dim(project_data_reduced) #2441x83 #remove percentage signs ir <- which(colnames(project_data_reduced)=="int_rate") for (row in 1:dim(project_data_reduced)[1]) { percent_sign <- regexpr("%", project_data_reduced[row,ir]) percent1 <- as.double(substr(project_data_reduced[row,ir],1,percent_sign-1)) project_data_reduced[row,ir] <- percent1 } project_data_reduced[,ir] <- sapply(project_data_reduced[,ir], as.double) ru <- which(colnames(project_data_reduced)=="revol_util") for (row in 1:dim(project_data_reduced)[1]) { percent_sign <- regexpr("%", project_data_reduced[row,ru]) percent2 <- as.double(substr(project_data_reduced[row,ru],1,percent_sign-1)) project_data_reduced[row,ru] <- percent2 } project_data_reduced[,ru] <- sapply(project_data_reduced[,ru], as.double) dim(project_data_reduced) #2441x83 #set "term" to be a single number ter <- which(colnames(project_data_reduced)=="term") for (row in 1:dim(project_data_reduced)[1]) { months_string_index <- regexpr('months', project_data_reduced[row,ter]) project_data_reduced[row,ter] <- substr(project_data_reduced[row,ter],2, months_string_index-2) } project_data_reduced[,ter] <- sapply(project_data_reduced[,ter], as.integer) dim(project_data_reduced) #2441x83 #Further removing of columns #These columns are linear combos of other columns #These should be deleted. del_cols_2 <- c("funded_amnt","funded_amnt_inv", "pymnt_plan", "recoveries", "application_type", "acc_now_delinq", "delinq_amnt", "num_tl_120dpd_2m", "num_tl_30dpd", "tax_liens", "hardship_flag", "debt_settlement_flag") project_data_reduced <- project_data_reduced[,! (names(project_data_reduced) %in% del_cols_2)] dim(project_data_reduced) #2441x71 project_data_reduced <- na.omit(project_data_reduced) dim(project_data_reduced) #1818x71  I started with a regular linear regression.  pd.lm.reg <- lm(formula = int_rate~., data = project_data_reduced)  And then did a semi-manual calculation of the MSE.  y.lm.reg <- predict(pd.lm.reg) sum((y.lm.reg - project_data_reduced$$int_rate)^2) #13580.01149 mean((y.lm.reg - project_data_reduced$$int_rate)^2) #7.469753293  The calculated mean square error is 7.469753293. Then, I wanted to show the usefulness of cross-validation, so I ran the leave-one-out-cross-validation (LOOCV) in R.  pd.glm.fit <- glm(int_rate~.,data=project_data_reduced) #This is the standard linear regression library(boot) set.seed(42) pd.cv.err.1818 <- cv.glm(data = project_data_reduced,pd.glm.fit, K=1818) pd.cv.err.1818 pd.cv.err.1818$delta


This gives the values

    #[1] 657,733,987.4 657,372,197.6


If I do a 10-fold validation like this

    set.seed(42)
pd.cv.err.10 <- cv.glm(data = project_data_reduced,pd.glm.fit, K=10)
pd.cv.err.10
pd.cv.err.10$delta  Then I get delta values of  #[1] 771,478,014.5 694,245,342.8  My understanding of the above code is that without specifying the cost function it is the case that the cv.glm function will automatically use an MSE calculation. In addition the documentation says that the \$delta term is the cross-validation estimate of the error. But the cross-validation estimates are roughly 8 orders of magnitude (10^8) times bigger than the calculated MSE. That doesn't seem right. Even if they were giving the total error rather than the mean then the mean values would be

    > pd.cv.err.1818$$delta/1818 [1] 361789.8721 361590.8677 > pd.cv.err.10$$delta/10
[1] 77147801.45 69424534.28


Why am I getting so vastly different numbers when doing the calculation manually vs. asking the cross-validation algorithm to look at it?

Am I misusing the code?

Am I calculating something incorrectly?

Does the $delta term actually mean something different than the MSE estimate after cross-validation? • Can you please either provide the code that was used to create your original sample or just give a new reproducible example? – usεr11852 says Reinstate Monic Oct 12 '18 at 23:10 • As requested I've added the 143 lines of cleaning code that I used. I'm not really interested in whether that code is clean or efficient or elegant or whatever other term you want to use. I just want to know whether I'm miscalculating or there's something wrong with the cv.glm command or what. – THill3 Oct 14 '18 at 20:55 • Thank you for this (+1) to your question. I do not get any crazy numbers like that from the $delta field; more like numbers ~7.7 which are indeed much closer to your "semi-manual calculation of the MSE". I suspect that something might be wrong with your R installation. Can you please reinstall R? – usεr11852 says Reinstate Monic Oct 14 '18 at 22:19
• Reinstalled R to version 3.5.1 and RStudio to version 1.1.456. That didn't fix it. Updated all of the packages using the RStudio functionality. That didn't fix it. Deleted the 'boot' package library (the actual files on the computer) and reinstalled it through RStudio (it used cran.rstudio.com/bin/windows/contrib/3.5/boot_1.3-20.zip). That didn't help. Tried running just the code in this post (perhaps some other code is messing with things?). That didn't help. Are you using the same version of the 'boot' package that I am using? – THill3 Oct 17 '18 at 12:03