Pred_on_BN_1236.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "With R kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Set the library path\n",
    "#.libPaths(\"/user/emma.foessing01/u11969/new_R_libs\")\n",
    "Sys.setenv(\"PKG_CXXFLAGS\"=\"-std=c++14\")\n",
    "\n",
    "print(R.version.string)\n",
    "\n",
    "# List of required packages\n",
    "list_of_packages <- c(\n",
    "  \"synthpop\", \"jsonlite\", \"codetools\", \"insight\", \"party\", \"haven\", \"dplyr\", \"rpart\", \"rpart.plot\",\n",
    "  \"randomForest\", \"pROC\", \"caret\", \"pracma\", \"here\", \"Hmisc\", \"purrr\",\n",
    "  \"ranger\", \"bnlearn\", \"arulesCBA\", \"network\", \"igraph\", \"xgboost\",\n",
    "  \"data.table\", \"doParallel\", \"parallel\", \"ExtDist\"\n",
    ")\n",
    "\n",
    "# Function to load packages and handle errors\n",
    "load_if_installed <- function(p) {\n",
    "  tryCatch({\n",
    "    library(p, character.only = TRUE)\n",
    "  }, error = function(e) {\n",
    "    message(sprintf(\"Package '%s' is not installed.\", p))\n",
    "  })\n",
    "}\n",
    "\n",
    "# Load all required packages\n",
    "lapply(list_of_packages, load_if_installed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "load(file = (paste0(here(), \"/cpspop.RData\")))\n",
    "cpspop <- cpspop[, c(setdiff(names(cpspop), c(\"income\", \"race\", \"marital\", \"educ\")), \"income\", \"race\", \"marital\", \"educ\")] #\n",
    "\n",
    "adult <- read.csv(file = (paste0(here(),\"/adult_preprocessed.csv\")))\n",
    "# delete NAs\n",
    "adult[adult == \"?\"] <- NA\n",
    "adult <- na.omit(adult)\n",
    "\n",
    "adult$workclass <- as.factor(adult$workclass)\n",
    "adult$education <- as.factor(adult$education)\n",
    "adult$marital_status <- as.factor(adult$marital_status)\n",
    "adult$relationship <- as.factor(adult$relationship)\n",
    "adult$race <- as.factor(adult$race)\n",
    "adult$sex <- as.factor(adult$sex)\n",
    "adult$native_country <- as.factor(adult$native_country)\n",
    "adult$income <- as.factor(adult$income)\n",
    "adult$occupation <- as.factor(adult$occupation)\n",
    "\n",
    "adult <- adult[, c(\"age\", \"fnlwgt\", \"capital_gain\", \"capital_loss\", \"hours_per_week\", \"income\", \"sex\", \"race\", \"relationship\", \"marital_status\", \"workclass\", \"occupation\", \"education\", \"native_country\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ein Datensatz zu genieren (m = 1) ist ausreichend, da ich keine Varianzanalyse machen werde. Damit die Ergebnisse nicht von einem zufälligen Prozess abhängen ist es sinnvoll über ein paar runs Mittelwerte zu bilden (50–100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "source": [
    "### Evaluation Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "## Calculate evaluation metrics for continuous targets\n",
    "evaluation_metrics_cont <- function(predictions, test_set){\n",
    "    # Residuals\n",
    "    residuals <- predictions - test_set$income\n",
    "    \n",
    "    # Mean Absolute Error (MAE)\n",
    "    MAE <- mean(abs(residuals))\n",
    "    \n",
    "    # Mean Squared Error (MSE) and Root Mean Squared Error (RMSE)\n",
    "    MSE <- mean(residuals^2)\n",
    "    RMSE <- sqrt(MSE)\n",
    "    \n",
    "    # R-squared: Guarding against zero variance in the target\n",
    "    SS_res <- sum(residuals^2)\n",
    "    SS_tot <- sum((test_set$income - mean(test_set$income))^2)\n",
    "    R_squared <- ifelse(SS_tot == 0, NA, 1 - (SS_res / SS_tot))\n",
    "    \n",
    "    # Mean Absolute Percentage Error (MAPE): Handling division by zero\n",
    "    MAPE <- ifelse(any(test_set$income == 0), NA, mean(abs(residuals / test_set$income)) * 100)\n",
    "    \n",
    "    metrics_df <- data.frame(\n",
    "        MAE = MAE, \n",
    "        MSE = MSE, \n",
    "        RMSE = RMSE, \n",
    "        R_squared = R_squared, \n",
    "        MAPE = MAPE\n",
    "    )\n",
    "    \n",
    "    return(metrics_df)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "## Calculate evaluation metrics for factored targets\n",
    "evaluation_metrics_factor <- function(predictions, test_set) {\n",
    "    # Ensure test_set is a data frame\n",
    "    test_set <- as.data.frame(test_set)\n",
    "    \n",
    "    # Ensure both predictions and test_set$income are factors with the same levels\n",
    "    predictions <- as.factor(predictions)\n",
    "    reference <- as.factor(test_set$income)\n",
    "    \n",
    "    # Ensure levels match between predictions and reference\n",
    "    levels(predictions) <- levels(reference)\n",
    "    \n",
    "    # Confusion matrix for the prediction on original data\n",
    "    cm <- caret::confusionMatrix(predictions, reference, mode = \"everything\")\n",
    "\n",
    "    # Saving evaluation metrics\n",
    "    accuracy <- cm$overall['Accuracy']\n",
    "    \n",
    "    if (length(levels(reference)) == 2) {\n",
    "        # Binary classification\n",
    "        f1 <- cm$byClass['F1']\n",
    "        sens <- cm$byClass['Sensitivity']\n",
    "        spec <- cm$byClass['Specificity']\n",
    "    } else {\n",
    "        # Multi-class classification: calculate metrics for each class and take the mean\n",
    "        f1 <- mean(cm$byClass[,'F1'], na.rm = TRUE)\n",
    "        sens <- mean(cm$byClass[,'Sensitivity'], na.rm = TRUE)\n",
    "        spec <- mean(cm$byClass[,'Specificity'], na.rm = TRUE)\n",
    "    }\n",
    "\n",
    "    # Create the dataframe\n",
    "    metrics_df <- data.frame(\n",
    "        Accuracy = accuracy, \n",
    "        F1 = f1, \n",
    "        Sensitivity = sens, \n",
    "        Specificity = spec\n",
    "    )\n",
    "    \n",
    "    return(metrics_df)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### CART"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cart_pred <- function(data, outer_folds, cp_steps, inner_folds) {\n",
    "    # Adjust evaluation metric to fit both numeric and factored targets\n",
    "    summaryFunctionType <- if (is.numeric(data$income)) defaultSummary else multiClassSummary\n",
    "\n",
    "    # Set control args for inner and outer loops\n",
    "    outer_control <- trainControl(method = \"cv\", number = outer_folds,\n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "    \n",
    "    inner_control <- trainControl(method = \"cv\", number = inner_folds, \n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "\n",
    "    # Define the grid for hyperparameter tuning\n",
    "    complexity <- 10^seq(log10(0.0001), log10(0.01), length.out = cp_steps)\n",
    "    tunegrid <- expand.grid(cp = complexity)\n",
    "\n",
    "    # Create outer CV folds\n",
    "    outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "\n",
    "    # Setup parallel backend for outer folds\n",
    "    num_cores <- detectCores() - 1  # Use all but one core\n",
    "    cl <- makeCluster(num_cores)\n",
    "    registerDoParallel(cl)\n",
    "\n",
    "    # Use foreach to parallelize the outer loop\n",
    "    outer_results <- foreach(i = seq_along(outer_cv_folds), .packages = c(\"caret\", \"dplyr\", \"rpart\"),\n",
    "                             .export = c(\"evaluation_metrics_cont\", \"evaluation_metrics_factor\")) %dopar% {\n",
    "        # Split data into outer folds\n",
    "        outer_test_index <- outer_cv_folds[[i]]\n",
    "        outer_testData <- data[outer_test_index, ]\n",
    "        outer_trainData <- data[-outer_test_index, ]\n",
    "\n",
    "        # Ensure the target variable in train and test sets have the same factor levels\n",
    "        outer_trainData$income <- factor(outer_trainData$income, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Add missing factor levels in income by temporarily adding rows with NA values for missing levels\n",
    "        missing_levels <- setdiff(levels(data$income), levels(outer_trainData$income))\n",
    "        if (length(missing_levels) > 0) {\n",
    "            for (level in missing_levels) {\n",
    "                # Add a temporary row with NA values and the missing level in income\n",
    "                outer_trainData <- rbind(outer_trainData, data.frame(income = factor(level, levels = levels(data$income))))\n",
    "            }\n",
    "        }\n",
    "\n",
    "        # Remove unused levels in all factor variables in training data\n",
    "        outer_trainData <- droplevels(outer_trainData)\n",
    "\n",
    "        # Hyperparameter tuning using inner CV\n",
    "        model <- caret::train(income ~ ., \n",
    "                              data = outer_trainData, \n",
    "                              method = \"rpart\", \n",
    "                              tuneGrid = tunegrid, \n",
    "                              trControl = inner_control,\n",
    "                              control = rpart.control(maxsurrogate = 0, maxcompete = 1))\n",
    "\n",
    "        # Get the best hyperparameters\n",
    "        best_hyperparameters <- model$bestTune\n",
    "\n",
    "        # Train the final model on the outer training set with the best hyperparameters\n",
    "        final_model <- caret::train(income ~ ., \n",
    "                                    data = outer_trainData, \n",
    "                                    method = \"rpart\", \n",
    "                                    trControl = outer_control, \n",
    "                                    tuneGrid = best_hyperparameters)\n",
    "\n",
    "        # Testing the final model on the outer test set\n",
    "        predictions <- predict(final_model, newdata = outer_testData)\n",
    "        \n",
    "        # Ensure consistent levels in predictions and test data for evaluation\n",
    "        predictions <- factor(predictions, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Evaluate the predictions\n",
    "        if (is.numeric(data$income)) {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData)\n",
    "        } else if (is.factor(data$income)) {\n",
    "            eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "        } else {\n",
    "            stop(\"The predicted target has to be numeric or factor.\")\n",
    "        }\n",
    "\n",
    "        return(eval)\n",
    "    }\n",
    "\n",
    "    stopCluster(cl)  # Stop the cluster when done\n",
    "\n",
    "    # Average the evaluation metrics over the outer folds\n",
    "    eval_avg_outer_folds <- do.call(rbind, outer_results) %>% as.data.frame() %>%\n",
    "                        dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "\n",
    "    return(eval_avg_outer_folds)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rf_pred <- function(data, outer_folds, mtry_steps, ntree_steps, inner_folds) {\n",
    "    # Adjust evaluation metric to fit both numeric and factored targets\n",
    "    summaryFunctionType <- if (is.numeric(data$income)) defaultSummary else multiClassSummary\n",
    "\n",
    "    # Set control args for inner and outer loops\n",
    "    outer_control <- trainControl(method = \"cv\", number = outer_folds,\n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "    \n",
    "    inner_control <- trainControl(method = \"cv\", number = inner_folds, \n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "\n",
    "    # Define the parameter grid for tuning\n",
    "    splitrule_value <- if (is.numeric(data$income)) \"variance\" else \"gini\"\n",
    "    \n",
    "    tunegrid <- expand.grid(mtry = seq(1, ncol(data) - 1, length.out = mtry_steps),\n",
    "                            splitrule = splitrule_value,\n",
    "                            min.node.size = if (is.numeric(data$income)) seq(5, 15, length.out = 10) else seq(1, 10, length.out = 10))\n",
    "\n",
    "    # Create outer CV folds\n",
    "    outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "\n",
    "    # Setup parallel backend for outer folds\n",
    "    num_cores <- detectCores() - 1  # Use all but one core\n",
    "    cl <- makeCluster(num_cores)\n",
    "    registerDoParallel(cl)\n",
    "\n",
    "    # Use foreach to parallelize the outer loop\n",
    "    outer_results <- foreach(i = seq_along(outer_cv_folds), .packages = c(\"caret\", \"dplyr\", \"ranger\"),\n",
    "                             .export = c(\"evaluation_metrics_cont\", \"evaluation_metrics_factor\")) %dopar% {\n",
    "        \n",
    "        # Split data into outer folds\n",
    "        outer_test_index <- outer_cv_folds[[i]]\n",
    "        outer_testData <- data[outer_test_index, ]\n",
    "        outer_trainData <- data[-outer_test_index, ]\n",
    "\n",
    "        # Ensure factor levels in the target variable for both training and test sets\n",
    "        outer_trainData$income <- factor(outer_trainData$income, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Add missing factor levels in income by temporarily adding rows with NA values for missing levels\n",
    "        missing_levels <- setdiff(levels(data$income), levels(outer_trainData$income))\n",
    "        if (length(missing_levels) > 0) {\n",
    "            for (level in missing_levels) {\n",
    "                # Add a temporary row with NA values and the missing level in income\n",
    "                outer_trainData <- rbind(outer_trainData, data.frame(income = factor(level, levels = levels(data$income))))\n",
    "            }\n",
    "        }\n",
    "\n",
    "        # Remove unused levels in all factor variables in training data\n",
    "        outer_trainData <- droplevels(outer_trainData)\n",
    "\n",
    "        # Hyperparameter tuning using inner CV\n",
    "        model <- caret::train(income ~ ., \n",
    "                              data = outer_trainData, \n",
    "                              method = \"ranger\",  \n",
    "                              tuneGrid = tunegrid, \n",
    "                              trControl = inner_control)\n",
    "\n",
    "        # Store the best hyperparameters\n",
    "        best_hyperparameters <- model$bestTune\n",
    "\n",
    "        # Train the final model on the outer training set with the best hyperparameters\n",
    "        final_model <- caret::train(income ~ ., \n",
    "                                    data = outer_trainData, \n",
    "                                    method = \"ranger\", \n",
    "                                    trControl = outer_control, \n",
    "                                    tuneGrid = best_hyperparameters)\n",
    "\n",
    "        # Test the final model on the outer test set\n",
    "        predictions <- predict(final_model, newdata = outer_testData)\n",
    "        \n",
    "        # Ensure consistent levels in predictions and test data for evaluation\n",
    "        predictions <- factor(predictions, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Evaluate the predictions and return a single-row data frame\n",
    "        if (is.numeric(data$income)) {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData)\n",
    "        } else if (is.factor(data$income)) {\n",
    "            eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "        } else {\n",
    "            stop(\"The predicted target has to be numeric or factor.\")\n",
    "        }\n",
    "\n",
    "        # Return eval as a single-row data frame\n",
    "        return(as.data.frame(t(unlist(eval))))\n",
    "    }\n",
    "\n",
    "    stopCluster(cl)  # Stop the cluster when done\n",
    "\n",
    "    # Combine outer_results into a data frame, then average\n",
    "    eval_df <- do.call(rbind, outer_results)\n",
    "\n",
    "    # Average the evaluation metrics over the outer folds\n",
    "    eval_avg_outer_folds <- eval_df %>%\n",
    "                            dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "\n",
    "    # Return the average evaluation metrics\n",
    "    return(eval_avg_outer_folds)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### XGB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "problematic datatype in xboost handling that was not compatible with caret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "xgb_pred <- function(data, outer_folds, inner_folds) {\n",
    "    # Save the original factor levels of the target variable (income) before conversion\n",
    "    original_levels <- levels(data$income)\n",
    "    \n",
    "    # Create outer folds once before the loop\n",
    "    outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "    \n",
    "    # Define the grid of hyperparameters to tune\n",
    "    tunegrid <- expand.grid(\n",
    "        max_depth = c(3, 5, 7),\n",
    "        eta = c(0.05, 0.1),\n",
    "        gamma = c(0, 0.1),\n",
    "        subsample = 0.8,\n",
    "        colsample_bytree = 0.8\n",
    "    )\n",
    "    \n",
    "    outer_results <- list()\n",
    "    \n",
    "    # Outer cross-validation loop using a regular for loop\n",
    "    for (i in seq_along(outer_cv_folds)) {\n",
    "        # Split data into outer folds\n",
    "        outer_test_index <- outer_cv_folds[[i]]\n",
    "        outer_testData <- data[outer_test_index, ]\n",
    "        outer_trainData <- data[-outer_test_index, ]\n",
    "    \n",
    "        # Align factor levels between training and test datasets\n",
    "        factor_vars <- names(Filter(is.factor, data))\n",
    "        \n",
    "        for (var in factor_vars) {\n",
    "            combined_levels <- union(levels(outer_trainData[[var]]), levels(outer_testData[[var]]))\n",
    "            outer_trainData[[var]] <- factor(outer_trainData[[var]], levels = combined_levels)\n",
    "            outer_testData[[var]]   <- factor(outer_testData[[var]], levels = combined_levels)\n",
    "        }\n",
    "        \n",
    "        # Remove unused factor levels in the training data\n",
    "        outer_trainData <- outer_trainData %>%\n",
    "            dplyr::mutate(across(where(is.factor), droplevels))\n",
    "\n",
    "        # Prepare the training and validation matrices\n",
    "        train_X <- data.matrix(outer_trainData[, !colnames(outer_trainData) %in% 'income'])\n",
    "        train_y <- outer_trainData$income\n",
    "        val_X <- data.matrix(outer_testData[, !colnames(outer_testData) %in% 'income'])\n",
    "        val_y <- outer_testData$income\n",
    "    \n",
    "        # Print debugging information\n",
    "        message(\"train_X dimensions: \", paste(dim(train_X), collapse = \" x \"))\n",
    "        message(\"train_y length: \", length(train_y))\n",
    "        message(\"val_X dimensions: \", paste(dim(val_X), collapse = \" x \"))\n",
    "        message(\"val_y length: \", length(val_y))\n",
    "    \n",
    "        if (is.factor(train_y)) {\n",
    "            # Convert factor levels to numeric starting from 0\n",
    "            train_y <- as.numeric(train_y) - 1\n",
    "            val_y <- as.numeric(val_y) - 1\n",
    "        }\n",
    "    \n",
    "        # Handle missing values\n",
    "        train_X[is.na(train_X)] <- 0\n",
    "        val_X[is.na(val_X)] <- 0\n",
    "    \n",
    "        # Create the DMatrix required for xgboost\n",
    "        dtrain <- xgb.DMatrix(data = train_X, label = train_y)\n",
    "        dtest <- xgb.DMatrix(data = val_X, label = val_y)\n",
    "    \n",
    "        best_params <- list()  # Initialize as an empty list\n",
    "        best_eval_metric <- Inf\n",
    "        best_nrounds <- 200  # Increased default value for nrounds\n",
    "    \n",
    "        # Iterate over all combinations of hyperparameters in tunegrid\n",
    "        for (j in 1:nrow(tunegrid)) {\n",
    "            # Determine the objective and eval_metric based on the target variable\n",
    "            if (is.numeric(data$income)) {\n",
    "                objective_name <- 'reg:squarederror'\n",
    "                eval_metric_name <- 'rmse'\n",
    "            } else if (length(unique(data$income)) == 2) {\n",
    "                objective_name <- 'binary:logistic'\n",
    "                eval_metric_name <- 'logloss'\n",
    "            } else {\n",
    "                objective_name <- 'multi:softprob'\n",
    "                eval_metric_name <- 'mlogloss'\n",
    "            }\n",
    "    \n",
    "            params <- list(\n",
    "                max_depth = tunegrid$max_depth[j],\n",
    "                eta = tunegrid$eta[j],\n",
    "                gamma = tunegrid$gamma[j],\n",
    "                subsample = tunegrid$subsample[j],\n",
    "                colsample_bytree = tunegrid$colsample_bytree[j],\n",
    "                objective = objective_name,\n",
    "                eval_metric = eval_metric_name\n",
    "            )\n",
    "    \n",
    "            # Include num_class only for multi-class classification\n",
    "            if (objective_name == 'multi:softprob') {\n",
    "                params$num_class <- length(unique(train_y))\n",
    "            }\n",
    "    \n",
    "            # Perform inner cross-validation using xgb.cv\n",
    "            cv_model <- xgb.cv(\n",
    "                params = params,\n",
    "                data = dtrain,\n",
    "                nrounds = 200,  # Increased number of boosting rounds\n",
    "                nfold = inner_folds,\n",
    "                early_stopping_rounds = 10,\n",
    "                verbose = FALSE\n",
    "            )\n",
    "    \n",
    "            # Retrieve the best metric from cross-validation\n",
    "            metric_column <- paste0('test_', eval_metric_name, '_mean')\n",
    "            current_eval_metric <- min(cv_model$evaluation_log[[metric_column]])\n",
    "            best_iteration <- which.min(cv_model$evaluation_log[[metric_column]])\n",
    "    \n",
    "            # If this model is better, update the best parameters\n",
    "            if (current_eval_metric < best_eval_metric) {\n",
    "                best_eval_metric <- current_eval_metric\n",
    "                best_params <- params\n",
    "                best_nrounds <- best_iteration\n",
    "            }\n",
    "        }\n",
    "    \n",
    "        # Train the final model on the outer training set with the best hyperparameters\n",
    "        final_model <- xgboost(\n",
    "            params = best_params,\n",
    "            data = dtrain,\n",
    "            nrounds = best_nrounds,\n",
    "            verbose = FALSE\n",
    "        )\n",
    "    \n",
    "        # Predict on the outer test set\n",
    "        predictions <- predict(final_model, dtest)\n",
    "    \n",
    "        if (is.factor(data$income) && length(unique(data$income)) > 2) {\n",
    "            # Multi-class classification\n",
    "            predicted_class_indices <- max.col(matrix(predictions, ncol = best_params$num_class, byrow = TRUE)) - 1\n",
    "            predictions <- factor(predicted_class_indices, levels = 0:(length(original_levels) - 1), labels = original_levels)\n",
    "        } else if (is.factor(data$income) && length(unique(data$income)) == 2) {\n",
    "            # Binary classification\n",
    "            predicted_classes <- ifelse(predictions >= 0.5, 1, 0)\n",
    "            predictions <- factor(predicted_classes, levels = 0:1, labels = original_levels)\n",
    "        }\n",
    "    \n",
    "        # Evaluate the model\n",
    "        if (is.factor(data$income)) {\n",
    "            eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "        } else {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData)\n",
    "        }\n",
    "    \n",
    "        # Store the evaluation result for this fold\n",
    "        outer_results[[i]] <- eval\n",
    "    }\n",
    "    \n",
    "    # Summarize the evaluation metrics across folds\n",
    "    eval_avg_outer_folds <- do.call(rbind, outer_results) %>%\n",
    "                            dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "    \n",
    "    return(eval_avg_outer_folds)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### BN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "this does not have an implemented function in the caret::train() function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "discretize_df = function(df, breaks = 5) {\n",
    "  for (var in colnames(df)) {\n",
    "    # Check if the variable is not a factor\n",
    "    if (!is.factor(df[[var]])) {\n",
    "\n",
    "      # Count the frequency of each unique value\n",
    "      freq_table <- table(df[[var]])\n",
    "\n",
    "      # Calculate the proportion of zeros, ensuring NA is handled\n",
    "      zero_proportion <- ifelse(!is.na(freq_table[as.character(0)]), \n",
    "                                freq_table[as.character(0)] / sum(freq_table), \n",
    "                                0)\n",
    "\n",
    "      # Determine the number of breaks based on zero proportion\n",
    "      if (zero_proportion > 4/5) {\n",
    "        new_breaks = 1\n",
    "      } else if (zero_proportion > 1/4) {\n",
    "        new_breaks = breaks - 2\n",
    "      } else if (zero_proportion > 1/5) {\n",
    "        new_breaks = breaks - 1\n",
    "      } else {\n",
    "        new_breaks = breaks\n",
    "      }\n",
    "      \n",
    "      # Separate zeros and non-zeros\n",
    "      zero_portion = (df[[var]] == 0)\n",
    "      non_zero_values = df[[var]][!zero_portion]\n",
    "\n",
    "      # Discretize non-zero values\n",
    "      if (length(non_zero_values) > 0) {\n",
    "        # Calculate breaks for non-zero values\n",
    "        range_values = range(non_zero_values, na.rm = TRUE)\n",
    "        breaks_values = seq(range_values[1], range_values[2], length.out = new_breaks + 1)\n",
    "        \n",
    "        # Ensure correct number of labels are created\n",
    "        labels = sapply(1:(length(breaks_values)-1), function(i) \n",
    "                        paste(\"(\", breaks_values[i], \"-\", breaks_values[i+1], \"]\", sep=\"\"))\n",
    "\n",
    "        # Use cut to apply these breaks and labels\n",
    "        discretized_non_zeros = cut(non_zero_values, breaks = breaks_values, labels = labels, include.lowest = TRUE)\n",
    "        # Combine zero and discretized non-zeros into the original dataframe\n",
    "        df[[var]] <- factor(ifelse(zero_portion, \"0\", as.character(discretized_non_zeros)))\n",
    "      } else {\n",
    "        # If all values are zero or the number of breaks is zero or negative\n",
    "        df[[var]] <- factor(\"0\")\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "  return(df)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "inner_cv <- function(data, target_var, folds, algorithms) {\n",
    "  # Create inner folds\n",
    "  inner_folds <- createFolds(data[[target_var]], k = folds)\n",
    "  \n",
    "  best_model <- NULL\n",
    "  best_performance <- -Inf\n",
    "  best_algorithm <- NULL\n",
    "  \n",
    "  for (algorithm in algorithms) {\n",
    "    cat(\"Trying algorithm:\", algorithm, \"\\n\")\n",
    "    fold_results <- c()\n",
    "    \n",
    "    for (i in seq_along(inner_folds)) {\n",
    "      inner_test_index <- inner_folds[[i]]\n",
    "      inner_trainData <- data[-inner_test_index, ]\n",
    "      inner_testData <- data[inner_test_index, ]\n",
    "      \n",
    "      # Fit Bayesian Network model using bnlearn algorithm\n",
    "      bn_model <- do.call(get(algorithm, envir = asNamespace(\"bnlearn\")), list(inner_trainData))\n",
    "      \n",
    "      # Fit the model to the training data\n",
    "      fitted_bn_model <- bnlearn::bn.fit(bn_model, inner_trainData)\n",
    "      \n",
    "      # Use Bayesian Likelihood Weighting for prediction\n",
    "      predictions <- predict(fitted_bn_model, node = target_var, data = inner_testData, method = \"bayes-lw\")\n",
    "      \n",
    "      # Handle missing levels in prediction\n",
    "      predictions <- factor(predictions, levels = levels(inner_trainData[[target_var]]))\n",
    "      \n",
    "      # Calculate the accuracy\n",
    "      accuracy <- mean(predictions == inner_testData[[target_var]], na.rm = TRUE)\n",
    "      fold_results[i] <- accuracy\n",
    "    }\n",
    "    \n",
    "    # Average performance for this algorithm\n",
    "    avg_performance <- mean(fold_results, na.rm = TRUE)\n",
    "    \n",
    "    if (!is.na(avg_performance) && avg_performance > best_performance) {\n",
    "      best_performance <- avg_performance\n",
    "      best_model <- fitted_bn_model\n",
    "      best_algorithm <- algorithm\n",
    "    }\n",
    "  }\n",
    "  \n",
    "  cat(\"Best algorithm selected:\", best_algorithm, \"with accuracy:\", best_performance, \"\\n\")\n",
    "  return(best_model)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "bn_pred <- function(data, outer_folds, inner_folds) {\n",
    "  # Discretize the data\n",
    "  data <- discretize_df(data)\n",
    "  \n",
    "  algorithms <- c(\"tabu\") # without hc\n",
    "  data$income <- factor(data$income, levels = unique(data$income))\n",
    "  \n",
    "  outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "  \n",
    "  # Setup parallel backend\n",
    "  num_cores <- detectCores() - 1  # You can adjust this depending on available cores\n",
    "  cl <- makeCluster(num_cores)\n",
    "  registerDoParallel(cl)\n",
    "  \n",
    "  # Use foreach for parallelization of outer folds\n",
    "  outer_results <- foreach(i = seq_along(outer_cv_folds), \n",
    "                           .packages = c(\"caret\", \"dplyr\", \"bnlearn\", \"doParallel\", \"foreach\", \"parallel\"), \n",
    "                           .export = c(\"discretize_df\", \"inner_cv\", \"evaluation_metrics_factor\")) %dopar% {\n",
    "    outer_test_index <- outer_cv_folds[[i]]\n",
    "    outer_testData <- data[outer_test_index, ]\n",
    "    outer_trainData <- data[-outer_test_index, ]\n",
    "    \n",
    "    # Get the best fitted model from inner CV\n",
    "    best_model <- inner_cv(outer_trainData, \"income\", inner_folds, algorithms)\n",
    "    \n",
    "    # Perform prediction using 'bayes-lw' method\n",
    "    predictions <- predict(best_model, node = \"income\", data = outer_testData, method = \"bayes-lw\")\n",
    "    \n",
    "    # Ensure both predictions and test data are factors\n",
    "    predictions <- factor(predictions, levels = levels(outer_testData$income))\n",
    "    \n",
    "    cat(\"Class of outer_testData: \", class(outer_testData), \"\\n\")\n",
    "\n",
    "    # Evaluate predictions (ensure test data is also a factor)\n",
    "    eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "    \n",
    "    return(eval)\n",
    "  }\n",
    "  \n",
    "  stopCluster(cl)  # Stop the cluster when done\n",
    "  registerDoSEQ()\n",
    "\n",
    "  # Average the evaluation metrics over the outer folds\n",
    "  eval_avg_outer_folds <- do.call(rbind, outer_results) %>%\n",
    "                          dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "  \n",
    "  return(eval_avg_outer_folds)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "svm_pred <- function(data, outer_folds, cost_steps, inner_folds) {\n",
    "    # Adjust evaluation metric to fit both numeric and factored targets\n",
    "    summaryFunctionType <- if (is.numeric(data$income)) defaultSummary else multiClassSummary\n",
    "\n",
    "    # Set control args for inner and outer loops\n",
    "    outer_control <- trainControl(method = \"cv\", number = outer_folds,\n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "    \n",
    "    inner_control <- trainControl(method = \"cv\", number = inner_folds, \n",
    "                                  summaryFunction = summaryFunctionType,\n",
    "                                  verboseIter = FALSE,\n",
    "                                  allowParallel = TRUE)\n",
    "\n",
    "    # Define the grid for hyperparameter tuning\n",
    "    cost_values <- 10^seq(log10(0.001), log10(100), length.out = cost_steps)\n",
    "    tunegrid <- expand.grid(C = cost_values, sigma = 0.1)  # sigma can also be tuned separately\n",
    "\n",
    "    # Create outer CV folds\n",
    "    outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "\n",
    "    # Setup parallel backend for outer folds\n",
    "    num_cores <- detectCores() - 1  # Use all but one core for parallel processing\n",
    "    cl <- makeCluster(num_cores)\n",
    "    registerDoParallel(cl)\n",
    "\n",
    "    # Parallelized outer loop: Cross-validation for model evaluation\n",
    "    outer_results <- foreach(i = seq_along(outer_cv_folds), .combine = rbind, .packages = c(\"caret\", \"dplyr\"),\n",
    "                             .export = c(\"evaluation_metrics_cont\", \"evaluation_metrics_factor\")) %dopar% {\n",
    "\n",
    "        # Split data into outer folds\n",
    "        outer_test_index <- outer_cv_folds[[i]]\n",
    "        outer_testData <- data[outer_test_index, ]\n",
    "        outer_trainData <- data[-outer_test_index, ]\n",
    "\n",
    "        # Ensure factor levels in the target variable for both training and test sets\n",
    "        outer_trainData$income <- factor(outer_trainData$income, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Add missing factor levels in income by temporarily adding rows with NA values for missing levels\n",
    "        missing_levels <- setdiff(levels(data$income), levels(outer_trainData$income))\n",
    "        if (length(missing_levels) > 0) {\n",
    "            for (level in missing_levels) {\n",
    "                # Add a temporary row with NA values and the missing level in income\n",
    "                outer_trainData <- rbind(outer_trainData, data.frame(income = factor(level, levels = levels(data$income))))\n",
    "            }\n",
    "        }\n",
    "\n",
    "        # Remove unused levels in all factor variables in training data\n",
    "        outer_trainData <- droplevels(outer_trainData)\n",
    "\n",
    "        # Hyperparameter tuning using inner CV\n",
    "        model <- caret::train(income ~ ., \n",
    "                              data = outer_trainData, \n",
    "                              method = \"svmRadial\",  \n",
    "                              tuneGrid = tunegrid, \n",
    "                              trControl = inner_control)\n",
    "\n",
    "        # Store the best hyperparameters\n",
    "        best_hyperparameters <- model$bestTune\n",
    "\n",
    "        # Train the final model on the outer training set with the best hyperparameters\n",
    "        final_model <- caret::train(income ~ ., \n",
    "                                    data = outer_trainData, \n",
    "                                    method = \"svmRadial\", \n",
    "                                    trControl = outer_control, \n",
    "                                    tuneGrid = best_hyperparameters)\n",
    "\n",
    "        # Testing the final model on the outer test set\n",
    "        predictions <- predict(final_model, newdata = outer_testData)\n",
    "\n",
    "        # Ensure consistent levels in predictions and test data for evaluation\n",
    "        predictions <- factor(predictions, levels = levels(data$income))\n",
    "        outer_testData$income <- factor(outer_testData$income, levels = levels(data$income))\n",
    "\n",
    "        # Evaluate the predictions and ensure it’s returned as a single-row data frame\n",
    "        if (is.numeric(data$income)) {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData)\n",
    "        } else if (is.factor(data$income)) {\n",
    "            eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "        } else {\n",
    "            stop(\"The predicted target has to be numeric or factor.\")\n",
    "        }\n",
    "\n",
    "        # Flatten the eval result to a consistent single-row data frame with explicit column names\n",
    "        eval <- as.data.frame(t(as.matrix(unlist(eval))))\n",
    "        colnames(eval) <- names(unlist(evaluation_metrics_factor(predictions, outer_testData)))\n",
    "\n",
    "        return(eval)\n",
    "    }\n",
    "\n",
    "    stopCluster(cl)  # Stop the cluster when done\n",
    "\n",
    "    # Combine outer_results into a data frame, then average\n",
    "    eval_avg_outer_folds <- outer_results %>%\n",
    "                            dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "\n",
    "    # Return the average evaluation metrics\n",
    "    return(eval_avg_outer_folds)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For BN synthezised data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "s <- 1236\n",
    "set.seed(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "simulation <- function(data, outer_folds = 5, inner_folds = 3,\n",
    "                       cp_steps = 10,  # CART params\n",
    "                       mtry_steps = 10, ntree_steps = 10,  # RF params\n",
    "                       nrounds_steps = 10, max_depth_steps = 10, eta_steps = 10, \n",
    "                       gamma_steps = 10, colsample_bytree_steps = 10, \n",
    "                       min_child_weight_steps = 10, subsample_steps = 10,  # XGB params\n",
    "                       cost_steps = 10  # SVM params\n",
    "                       ) {\n",
    "  \n",
    "  # Check if the target variable is numeric (regression) or factor (classification)\n",
    "  target_is_numeric <- is.numeric(data$income)\n",
    "  \n",
    "  # Initialize empty lists to store evaluation dataframes\n",
    "  eval_list <- list()  # For regression or classification metrics\n",
    "  eval_bn_list <- list()  # For BN metrics when the target is numeric (regression)\n",
    "  \n",
    "  # create synthetic data\n",
    "  syndata <- readRDS(paste0(here(), \"/results/\", deparse(substitute(data)), \"_bn_\", as.character(s), \".rds\"))\n",
    "  \n",
    "  # Prediction models with nested CV and grid search\n",
    "  cat(\"SVM eval\")\n",
    "  SVM_eval <- svm_pred(syndata, outer_folds, cost_steps, inner_folds)\n",
    "  cat(\"BN eval\")\n",
    "  BN_eval <- bn_pred(syndata, outer_folds, inner_folds)\n",
    "  cat(\"RF eval\")\n",
    "  RF_eval <- rf_pred(syndata, outer_folds, mtry_steps, ntree_steps, inner_folds)\n",
    "  cat(\"XGB eval\")\n",
    "  XGB_eval <- xgb_pred(syndata, outer_folds, inner_folds)\n",
    "  cat(\"CART eval\")\n",
    "  CART_eval <- cart_pred(syndata, outer_folds, cp_steps, inner_folds)\n",
    "  \n",
    "\n",
    "  cat(\"CART Evaluation:\\n\")\n",
    "  cat(capture.output(CART_eval), sep = \"\\n\")\n",
    "\n",
    "  cat(\"RF Evaluation:\\n\")\n",
    "  cat(capture.output(RF_eval), sep = \"\\n\")\n",
    "\n",
    "  cat(\"BN Evaluation:\\n\")\n",
    "  cat(capture.output(BN_eval), sep = \"\\n\")\n",
    "\n",
    "  cat(\"XGB Evaluation:\\n\")\n",
    "  cat(capture.output(XGB_eval), sep = \"\\n\")\n",
    "\n",
    "  cat(\"SVM Evaluation:\\n\")\n",
    "  cat(capture.output(SVM_eval), sep = \"\\n\")\n",
    "\n",
    "  rm(syndata)\n",
    "  gc()\n",
    "  \n",
    "  if (target_is_numeric) {\n",
    "    # Combine only regression metrics (CART, RF, XGB, SVM)\n",
    "    eval <- rbind(CART = CART_eval, RF = RF_eval, XGB = XGB_eval, SVM = SVM_eval)\n",
    "    #eval_list[[i]] <- eval\n",
    "    \n",
    "    # Save BN metrics separately for regression task\n",
    "    eval_bn <- BN_eval\n",
    "    #eval_bn_list[[i]] <- BN_eval\n",
    "\n",
    "    results <- list(eval, eval_bn)\n",
    "  } else {\n",
    "    # Combine all classification metrics, including BN\n",
    "    eval <- rbind(CART = CART_eval, RF = RF_eval, XGB = XGB_eval, SVM = SVM_eval, BN = BN_eval)\n",
    "    #eval_list[[i]] <- eval\n",
    "\n",
    "    results <- list(eval)\n",
    "  }\n",
    "  \n",
    "  return(results)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cps <- cpspop\n",
    "rm(cpspop)\n",
    "cps_res <- simulation(data = cps, outer_folds = 5, inner_folds = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "save(cps_res, file = paste0(here(), \"/results/cps_bn_res_\", as.character(s) ,\".RData\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#adult_res <- simulation(data = adult, outer_folds = 5, inner_folds = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#save(adult_res, file = paste0(here(), \"/results/adult_bn_res_\", as.character(s) ,\".RData\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#adult_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cps_res"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}