Thank you very much! Yes, I checked my code and found no fewer than three bugs. For the record:
- In
extract_log_lik_K_memory_friendly(), I should have the lines:
for(k in 1:K){
load(paste(save_file_stem,"k=",k,sep="")) # "merged_model"
instead of
for(k in 1:K){
load(paste(save_file_stem,"k=1",sep="")) # "merged_model"
That alone would produce the behavior you described, because 90% of the log likelihoods that I considered “holdout data” were actually computed from training data. Secondly, it turns out that in do_K_fold_cross_validation_memory_friendly(), I extracted the holdout likelihood correctly for the purpose of computing the expected log posterior densities, but (perhaps ironically) extracted it incorrectly when my goal was to return the data in a different format in order to pass into my get_loo_from_kfold() function for the purpose of catching bugs. Here’s the corrected code in case anyone else wants to use this:
do_K_fold_cross_validation_memory_friendly <- function(data,
stan_file,
save_file_stem,
do_stratification=TRUE,
n_folds=10,
n_chains=3,
n_cores=3,
#merge_chains=FALSE,
iter=2000,
control=list(max_treedepth = 15),
desired_stratification_percentage=0.015,
truncation_percentage=0.02,
debug=TRUE,
recalculate=FALSE)
{
library(rstan)
rstan_options(auto_write = TRUE)
library(pbmcapply)
library(loo)
library(caret)
K <- n_folds
if(!(recalculate)){
## TODO: Make this its own function. (Get holdout data in flattened and unflattened formats)
flattened.data <- as.vector(as.matrix(data$Y))
## For now, the stratification is done locally.
if(do_stratification){
ones_and_twos <- get_group_membership(data,
desired_group_1_percentage=desired_stratification_percentage,
truncation_percentage=truncation_percentage,
debug=TRUE)
}else{
ones_and_twos <- rep(1, length(flattened.data))
ones_and_twos[((length(flattened.data)/2)+1):length(flattened.data)] <- 2
}
folds <- createFolds(ones_and_twos, n_folds, list=FALSE)
# TODO: Testing my work: I want the number of "yes" responses in each fold to be about the same.
holdout_K_from_caret <- list()
for(k in 1:K)
{
holdout_vector <- rep(0,nrow(data$Y)*ncol(data$Y))
holdout_vector[which(folds==k)] <- 1
holdout_K_from_caret[[k]] <- unflatten(holdout_vector,nrow=nrow(data$Y),ncol=ncol(data$Y))
}
holdout_K_from_caret_flattened <- list()
for(i in 1:K){
holdout_K_from_caret_flattened[[i]] <- flatten(holdout_K_from_caret[[i]])
}
save(holdout_K_from_caret_flattened, file=paste(save_file_stem,"holdout_indices",sep=""))
#run the functions
I <- dim(data)[1]
J <- dim(data)[2]
data_m <- data
#create a list of data list
data_l <- rep(list(data),K)
#add the holdout index to it
for(i in 1:K) data_l[[i]]$holdout <- holdout_K_from_caret[[i]]
cat("In do_K_fold_cross_validation: Running K memory-friendly Stan models...\n")
stan_kfold_memory_friendly(file=stan_file,
save_file_stem=save_file_stem,
data_l,
chains=n_chains,
cores=n_cores,
iter=iter,
control=control)
# This doesn't return a list of Stan fit objects. So I want to pass in the
# save file stem to extract_log_lik_K.
}
if(recalculate){
load(paste(save_file_stem,"holdout_indices",sep="")) # holdout_K_from_caret_flattened
}
cat("In do_K_fold_cross_validation: extracting log-likelihood...\n")
ees <- extract_log_lik_K_memory_friendly(save_file_stem=save_file_stem,
n_folds=K,
list_of_holdout=holdout_K_from_caret_flattened)
ee_for_loo <- ees$log_lik_loo
ee_for_kfold <- ees$log_lik_kfold
#ee_for_loo <- extract_log_lik_K_memory_friendly(save_file_stem=save_file_stem,
# n_folds=K,
# list_of_holdout=holdout_K_from_caret_flattened,
# merge_chains=merge_chains,
# debug=debug)
#ee_for_kfold <- extract_log_lik_K_memory_friendly(save_file_stem=save_file_stem,
# n_folds=K,
# list_of_holdout=holdout_K_from_caret_flattened,
# debug=debug)
#r_eff <- relative_eff(exp(ee))
cat("In do_K_fold_cross_validation: calculating ELPD and standard error...\n")
kk <- kfold(ee_for_kfold)
result <- list()
result$kk <- kk
result$ee_for_loo <- ee_for_loo
result$ee_for_kfold <- ee_for_kfold
return(result)
}
#ee_for_loo <- extract_log_lik_K_memory_friendly(save_file_stem=save_file_stem,
# n_folds=K,
# list_of_holdout=holdout_K_from_caret_flattened,
# merge_chains=merge_chains,
# debug=debug)
#ee_for_kfold <- extract_log_lik_K_memory_friendly(save_file_stem=save_file_stem,
# n_folds=K,
# list_of_holdout=holdout_K_from_caret_flattened,
# debug=debug)
- Perhaps most importantly, even with the previous bugs fixed, my code still would not work because I have label switching taking place between chains. I’ve seen that before and addressed the problem, but with K-fold cross-validation I have the added responsibility of making sure that the labels I am attaching to each chain are consistent across the model runs. Because this was not the case in my previous analysis, it’s no wonder the 3-D model appears to have the best fit, because I’m essentially just fitting noise!
I’m working on some plotting functions now to make sure that I can better visualize the assignment of labels to chains across all ten folds, in order to make my analysis less error-prone. I’ll post back once I’ve got the code working, it may be a couple of weeks as I’ll have to put this project down to do some other work. Thank you very much for your assistance, Dr. Vehtari! The feedback definitely helped me to catch these bugs.
Gratefully,
Richard