Pred_on_orig.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "With R kernel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Further info\n",
    "Models usable with train() from caret: <br>\n",
    "https://topepo.github.io/caret/train-models-by-tag.html#Model_Tree <br>"
   ]
  },
  {
   "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\"\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": [
    "## Parallelisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Detect the number of available cores\n",
    "num_cores <- detectCores() - 1  # Leave one core free for system operations\n",
    "\n",
    "# Register the parallel backend\n",
    "cl <- makeCluster(num_cores)\n",
    "registerDoParallel(cl)"
   ]
  },
  {
   "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": [
    "# just the prediction\n",
    "cart_pred <- function(data, outer_folds, cp_steps, inner_folds){############adjust##############\n",
    "    # adjust evaluation metric to fit both numeric and factored targets\n",
    "    summaryFunctionType <- if (is.numeric(data$income)) defaultSummary else multiClassSummary\n",
    "    # metric: train() uses per default RSME and Accuracy for numeric and factored targets\n",
    "\n",
    "    #  set control args\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)############adjust##############\n",
    "\n",
    "    # Create grid\n",
    "    tunegrid <- expand.grid(cp = complexity)############adjust##############\n",
    "\n",
    "    # Initialize variables to store results\n",
    "    outer_results <- list()\n",
    "\n",
    "    outer_cv_folds = createFolds(data$income, k = outer_folds)\n",
    "    \n",
    "    # Outer loop: Cross-validation for model evaluation\n",
    "    for (i in seq_along(outer_folds)) {\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",
    "        # Hyperparameter tuning using inner CV\n",
    "        # No need for inner loop because \"train\" does k-fold CV already\n",
    "        model <- caret::train(income ~ ., \n",
    "                        data = outer_trainData, \n",
    "                        method = \"rpart\", ############adjust##############\n",
    "                        tuneGrid = tunegrid, \n",
    "                        trControl = inner_control,\n",
    "                        control = rpart.control(maxsurrogate = 0, maxcompete = 1) ############adjust##############\n",
    "                        )#,\n",
    "                        #metric = metricType)\n",
    "            \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 = \"rpart\",############adjust##############\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",
    "        if (is.numeric(data$income)) {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData) # postResample is a useful caret function\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",
    "        # Store the evaluation metrics for this outer fold\n",
    "        outer_results[[i]] <- eval\n",
    "    }\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",
    "    \n",
    "\n",
    "    # Return the average evaluation metrics\n",
    "    return(eval_avg_outer_folds)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### RF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Prediction function for Random Forest (Parallelized)\n",
    "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\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",
    "    # Initialize variables to store results\n",
    "    outer_cv_folds <- createFolds(data$income, k = outer_folds)\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",
    "        # 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",
    "        # Testing the final model on the outer test set\n",
    "        predictions <- predict(final_model, newdata = outer_testData)\n",
    "        \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 the evaluation metrics for this outer fold\n",
    "        return(eval)\n",
    "    }\n",
    "\n",
    "    # Average the evaluation metrics over the outer folds\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": [
    "##### 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",
    "\n",
    "    # Save the original factor levels of the target variable (income) before conversion\n",
    "    original_levels <- levels(data$income)\n",
    "\n",
    "    # Initialize variables to store results\n",
    "    outer_results <- list()\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",
    "    for (i in seq_along(outer_cv_folds)) {\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",
    "        # 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",
    "        if (is.factor(train_y)) {\n",
    "            # Ensure consistent factor levels\n",
    "            levels(train_y) <- levels(data$income)\n",
    "            levels(val_y) <- levels(data$income)\n",
    "            \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",
    "            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 = if (is.numeric(data$income)) 'reg:squarederror' else if (length(unique(data$income)) == 2) 'binary:logistic' else 'multi:softprob',\n",
    "                eval_metric = if (is.numeric(data$income)) 'rmse' else if (length(unique(data$income)) == 2) 'logloss' else 'mlogloss'\n",
    "            )\n",
    "\n",
    "            # Include num_class only for multi-class classification\n",
    "            if (is.factor(data$income) && length(unique(data$income)) > 2) {\n",
    "                params$num_class <- length(unique(train_y))\n",
    "            }\n",
    "\n",
    "            # Perform inner cross-validation using xgb.cv\n",
    "            message(\"Trying params: \", params)  # Log current params\n",
    "            cv_model <- xgb.cv(\n",
    "                params = params,\n",
    "                data = dtrain,\n",
    "                nrounds = 200,  # Increased number of boosting rounds\n",
    "                nfold = inner_folds,\n",
    "                verbose = FALSE\n",
    "            )\n",
    "\n",
    "            # Check if the cross-validation completed successfully\n",
    "            if (is.null(cv_model)) {\n",
    "                message(\"Inner CV failed for params: \", params)\n",
    "                next\n",
    "            }\n",
    "\n",
    "            # Retrieve the best metric from cross-validation\n",
    "            current_eval_metric <- min(cv_model$evaluation_log$test_rmse_mean)\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 <- which.min(cv_model$evaluation_log$test_rmse_mean)\n",
    "            }\n",
    "        }\n",
    "\n",
    "        # Ensure best_params is not empty\n",
    "        if (length(best_params) == 0) {\n",
    "            message(\"No valid parameters found, using default parameters.\")\n",
    "            best_params <- list(\n",
    "                max_depth = 3,\n",
    "                eta = 0.1,\n",
    "                gamma = 0,\n",
    "                subsample = 0.8,\n",
    "                colsample_bytree = 0.8,\n",
    "                objective = 'binary:logistic',\n",
    "                eval_metric = 'logloss'\n",
    "            )\n",
    "            best_nrounds <- 200  # Default value\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 = length(unique(train_y)), 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",
    "        } else if (is.numeric(data$income)) {\n",
    "            eval <- evaluation_metrics_cont(predictions, outer_testData)\n",
    "        } else {\n",
    "            stop(\"The predicted target has to be numeric or factor.\")\n",
    "        }\n",
    "\n",
    "        # Ensure predictions and true labels have the same length\n",
    "        if (length(predictions) != nrow(outer_testData)) {\n",
    "            stop(\"Length of predictions does not match number of test samples.\")\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",
    "        outer_results[[i]] <- eval\n",
    "    }\n",
    "\n",
    "    # Filter out NULL values from outer_results\n",
    "    valid_results <- Filter(Negate(is.null), outer_results)\n",
    "\n",
    "    # Check if there are valid results to summarize\n",
    "    if (length(valid_results) == 0) {\n",
    "        stop(\"No valid model results were obtained.\")\n",
    "    }\n",
    "\n",
    "    # Summarize the evaluation metrics across folds\n",
    "    eval_avg_outer_folds <- do.call(rbind, valid_results) %>%\n",
    "                            dplyr::summarise(across(everything(), mean, na.rm = TRUE))\n",
    "\n",
    "    return(eval_avg_outer_folds)\n",
    "\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(\"hc\", \"tabu\")  # You can add more algorithms here\n",
    "  data$income <- factor(data$income, levels = unique(data$income))\n",
    "\n",
    "  outer_results <- list()\n",
    "  outer_cv_folds <- createFolds(data$income, k = outer_folds)\n",
    "  \n",
    "  for (i in seq_along(outer_cv_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",
    "    # 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",
    "    eval <- evaluation_metrics_factor(predictions, outer_testData)\n",
    "    outer_results[[i]] <- eval\n",
    "  }\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",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Prediction function for SVM with Radial Kernel (Parallelized)\n",
    "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\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)  # Adjust the range as needed\n",
    "    tunegrid <- expand.grid(C = cost_values, sigma = 0.1)  # sigma can also be tuned separately\n",
    "\n",
    "    # Initialize variables to store results\n",
    "    outer_results <- list()\n",
    "    outer_cv_folds = createFolds(data$income, k = outer_folds)\n",
    "    \n",
    "    # Outer loop: Cross-validation for model evaluation (Parallelized)\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",
    "        # 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",
    "\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",
    "        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 the evaluation metrics for this outer fold\n",
    "        return(eval)\n",
    "    }\n",
    "\n",
    "    # Average the evaluation metrics over the outer folds\n",
    "    eval_avg_outer_folds <- outer_results %>%\n",
    "                            summarise(across(everything(), mean, na.rm = TRUE))\n",
    "\n",
    "    # Return the average evaluation metrics\n",
    "    return(eval_avg_outer_folds)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For original data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "simulation <- function(data, nrun = 10, 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 = 2  # 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",
    "  # Set initial seed\n",
    "  s <- 1234\n",
    "  for (i in 1:nrun) {\n",
    "    # Vary seed with each run\n",
    "    s <- s + 1\n",
    "    set.seed(s)\n",
    "    \n",
    "    # Prediction models with nested CV and grid search\n",
    "    cat(\"CART eval\")\n",
    "    CART_eval <- cart_pred(syndata, outer_folds, cp_steps, 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(\"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",
    "    \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_list[[i]] <- BN_eval\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",
    "    \n",
    "    cat(\"Run\", i, \"completed\")\n",
    "  }\n",
    "  \n",
    "  # Now average the results over all runs\n",
    "  if (target_is_numeric) {\n",
    "    # Average regression metrics\n",
    "    sum_regression_df <- Reduce(\"+\", eval_list)\n",
    "    eval_regression_avg <- sum_regression_df / nrun\n",
    "    \n",
    "    # Average BN classification metrics (separately)\n",
    "    sum_bn_df <- Reduce(\"+\", eval_bn_list)\n",
    "    eval_bn_avg <- sum_bn_df / nrun\n",
    "    \n",
    "    # Return both regression and BN classification metrics\n",
    "    results <- list(eval_regression_avg = eval_regression_avg, eval_bn_avg = eval_bn_avg)\n",
    "  } else {\n",
    "    # Average classification metrics (including BN)\n",
    "    sum_classification_df <- Reduce(\"+\", eval_list)\n",
    "    eval_classification_avg <- sum_classification_df / nrun\n",
    "    \n",
    "    # Return the classification average (including BN)\n",
    "    results <- list(eval_classification_avg = eval_classification_avg)\n",
    "  }\n",
    "  \n",
    "  return(results)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "rownames <- row.names(eval_bn_avg[[1]])\n",
    "        row.names(eval_bn_avg) <- rownames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "adult_res <- simulation(data = adult, nrun = 10, outer_folds = 5, inner_folds = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cps_res <- simulation(data = cpspop, nrun = 10, outer_folds = 5, inner_folds = 3)"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Stop the parallel cluster\n",
    "stopCluster(cl)\n",
    "registerDoSEQ()  # Return to sequential computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Save the data as .csv because RData files could not be opened by other R version\n",
    "save(cps_res, file = paste0(here(), \"/results/cps_orig_res.RData\"))\n",
    "save(adult_res, file = paste0(here(), \"/results/adult_orig_res.RData\"))"
   ]
  }
 ],
 "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
}