|
| 1 | +# |
| 2 | +# Copyright (C) 2013-2021 University of Amsterdam |
| 3 | +# |
| 4 | +# This program is free software: you can redistribute it and/or modify |
| 5 | +# it under the terms of the GNU General Public License as published by |
| 6 | +# the Free Software Foundation, either version 2 of the License, or |
| 7 | +# (at your option) any later version. |
| 8 | +# |
| 9 | +# This program is distributed in the hope that it will be useful, |
| 10 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 | +# GNU General Public License for more details. |
| 13 | +# |
| 14 | +# You should have received a copy of the GNU General Public License |
| 15 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. |
| 16 | +# |
| 17 | + |
| 18 | +mlClassificationLogistic <- function(jaspResults, dataset, options, ...) { |
| 19 | + |
| 20 | + # Preparatory work |
| 21 | + dataset <- .mlClassificationReadData(dataset, options) |
| 22 | + .mlClassificationErrorHandling(dataset, options, type = "logistic") |
| 23 | + |
| 24 | + # Check if analysis is ready to run |
| 25 | + ready <- .mlClassificationReady(options, type = "logistic") |
| 26 | + |
| 27 | + # Compute results and create the model summary table |
| 28 | + .mlClassificationTableSummary(dataset, options, jaspResults, ready, position = 1, type = "logistic") |
| 29 | + |
| 30 | + # If the user wants to add the classes to the data set |
| 31 | + .mlClassificationAddPredictionsToData(dataset, options, jaspResults, ready) |
| 32 | + |
| 33 | + # Add test set indicator to data |
| 34 | + .mlAddTestIndicatorToData(options, jaspResults, ready, purpose = "classification") |
| 35 | + |
| 36 | + # Create the data split plot |
| 37 | + .mlPlotDataSplit(dataset, options, jaspResults, ready, position = 2, purpose = "classification", type = "logistic") |
| 38 | + |
| 39 | + # Create the confusion table |
| 40 | + .mlClassificationTableConfusion(dataset, options, jaspResults, ready, position = 3) |
| 41 | + |
| 42 | + # Create the class proportions table |
| 43 | + .mlClassificationTableProportions(dataset, options, jaspResults, ready, position = 4) |
| 44 | + |
| 45 | + # Create the validation measures table |
| 46 | + .mlClassificationTableMetrics(dataset, options, jaspResults, ready, position = 5) |
| 47 | + |
| 48 | +# # Create the variable importance table |
| 49 | +# .mlTableFeatureImportance(options, jaspResults, ready, position = 6, purpose = "classification") |
| 50 | + |
| 51 | +# # Create the shap table |
| 52 | +# .mlTableShap(dataset, options, jaspResults, ready, position = 7, purpose = "classification") |
| 53 | + |
| 54 | +# # Create the ROC curve |
| 55 | +# .mlClassificationPlotRoc(dataset, options, jaspResults, ready, position = 8, type = "logistic") |
| 56 | + |
| 57 | + # Create the Andrews curves |
| 58 | + .mlClassificationPlotAndrews(dataset, options, jaspResults, ready, position = 9) |
| 59 | + |
| 60 | +# # Decision boundaries |
| 61 | +# .mlClassificationPlotBoundaries(dataset, options, jaspResults, ready, position = 10, type = "logistic") |
| 62 | +} |
| 63 | + |
| 64 | +.logisticRegressionClassification <- function(dataset, options, jaspResults, ready) { |
| 65 | + # Import model formula from jaspResults |
| 66 | + formula <- jaspResults[["formula"]]$object |
| 67 | + # Split the data into training and test sets |
| 68 | + if (options[["holdoutData"]] == "testSetIndicator" && options[["testSetIndicatorVariable"]] != "") { |
| 69 | + # Select observations according to a user-specified indicator (included when indicator = 1) |
| 70 | + trainingIndex <- which(dataset[, options[["testSetIndicatorVariable"]]] == 0) |
| 71 | + } else { |
| 72 | + # Sample a percentage of the total data set |
| 73 | + trainingIndex <- sample.int(nrow(dataset), size = ceiling((1 - options[["testDataManual"]]) * nrow(dataset))) |
| 74 | + } |
| 75 | + trainingSet <- dataset[trainingIndex, ] |
| 76 | + # Create the generated test set indicator |
| 77 | + testIndicatorColumn <- rep(1, nrow(dataset)) |
| 78 | + testIndicatorColumn[trainingIndex] <- 0 |
| 79 | + # Just create a train and a test set (no optimization) |
| 80 | + testSet <- dataset[-trainingIndex, ] |
| 81 | + if (nlevels(trainingSet[[options[["target"]]]]) == 2) { |
| 82 | + family = "binomial" |
| 83 | + trainingFit <- stats::glm(formula, data = trainingSet, family = family) |
| 84 | + # Use the specified model to make predictions for dataset |
| 85 | + testPredictions <- levels(trainingSet[[options[["target"]]]])[round(predict(trainingFit, newdata = testSet, type = "response"), 0) + 1] |
| 86 | + dataPredictions <- levels(trainingSet[[options[["target"]]]])[round(predict(trainingFit, newdata = dataset, type = "response"), 0) + 1] |
| 87 | + } else { |
| 88 | + family <- "multinomial" |
| 89 | + trainingFit <- VGAM::vglm(formula, data = trainingSet, family = family) |
| 90 | + # Use the specified model to make predictions for dataset |
| 91 | + testPredictions <- .mlClassificationMultinomialPredictions(trainingSet, options, predict(trainingFit, newdata = testSet)) |
| 92 | + dataPredictions <- .mlClassificationMultinomialPredictions(trainingSet, options, predict(trainingFit, newdata = dataset)) |
| 93 | + } |
| 94 | + # Create results object |
| 95 | + result <- list() |
| 96 | + result[["formula"]] <- formula |
| 97 | + result[["family"]] <- family |
| 98 | + result[["model"]] <- trainingFit |
| 99 | + result[["confTable"]] <- table("Pred" = testPredictions, "Real" = testSet[, options[["target"]]]) |
| 100 | + result[["testAcc"]] <- sum(diag(prop.table(result[["confTable"]]))) |
| 101 | +# result[["auc"]] <- .classificationCalcAUC(testSet, trainingSet, options, "logisticClassification") |
| 102 | + result[["ntrain"]] <- nrow(trainingSet) |
| 103 | + result[["ntest"]] <- nrow(testSet) |
| 104 | + result[["testReal"]] <- testSet[, options[["target"]]] |
| 105 | + result[["testPred"]] <- testPredictions |
| 106 | + result[["train"]] <- trainingSet |
| 107 | + result[["test"]] <- testSet |
| 108 | + result[["testIndicatorColumn"]] <- testIndicatorColumn |
| 109 | + result[["classes"]] <- dataPredictions |
| 110 | +# result[["explainer"]] <- DALEX::explain(result[["model"]], type = "classification", data = result[["train"]], y = result[["train"]][, options[["target"]]], predict_function = function(model, data) predict(model, newdata = data, type = "raw")) |
| 111 | +# if (nlevels(result[["testReal"]]) == 2) { |
| 112 | +# result[["explainer_fi"]] <- DALEX::explain(result[["model"]], type = "classification", data = result[["train"]], y = as.numeric(result[["train"]][, options[["target"]]]) - 1, predict_function = function(model, data) predict(model, newdata = data, type = "class")) |
| 113 | +# } else { |
| 114 | +# result[["explainer_fi"]] <- DALEX::explain(result[["model"]], type = "multiclass", data = result[["train"]], y = result[["train"]][, options[["target"]]] , predict_function = function(model, data) predict(model, newdata = data, type = "raw")) |
| 115 | +# } |
| 116 | + return(result) |
| 117 | +} |
| 118 | + |
| 119 | +.mlClassificationMultinomialPredictions <- function(trainingSet, options, predictions) { |
| 120 | + num_categories <- ncol(predictions) + 1 |
| 121 | + probs <- matrix(0, nrow = nrow(predictions), ncol = num_categories) |
| 122 | + for (i in 1:(num_categories - 1)) { |
| 123 | + probs[, i] <- exp(predictions[, i]) |
| 124 | + } |
| 125 | + probs[, num_categories] <- 1 |
| 126 | + row_sums <- rowSums(probs) |
| 127 | + probs <- probs / row_sums |
| 128 | + predicted_category <- apply(probs, 1, which.max) |
| 129 | + categories <- levels(trainingSet[[options[["target"]]]]) |
| 130 | + predicted_categories <- categories[predicted_category] |
| 131 | + return(predicted_categories) |
| 132 | +} |
0 commit comments