## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( dpi=300, fig.align = "center", out.width = "80%", echo = TRUE, message = FALSE, warning = FALSE ) ## ----message=FALSE, warning=FALSE--------------------------------------------- # install.packages("GeDS") library("GeDS") ## ----------------------------------------------------------------------------- # Generate a data sample for the response variable # Y and the single covariate X set.seed(123) N <- 500 f_1 <- function(x) (10*x/(1+100*x^2))*4+4 X <- sort(runif(N, min = -2, max = 2)) # Specify a model for the mean of Y to include only a component # non-linear in X, defined by the function f_1 means <- f_1(X) # Add (Normal) noise to the mean of Y Y <- rnorm(N, means, sd = 0.2) data = data.frame(X, Y) Gmodboost <- NGeDSboost(Y ~ f(X), data = data) ## ----fig.width=16, fig.height=8----------------------------------------------- layout(matrix(c(1,2), nrow=1, byrow=TRUE)) visualize_boosting(Gmodboost, 0:N.boost.iter(Gmodboost), final_fits = TRUE) ## ----------------------------------------------------------------------------- # Load the processed version of the Ames dataset library(AmesHousing) ames_data <- make_ames() # Keep only numeric variables ames_data <- ames_data[sapply(ames_data, is.numeric)] # Define the response variable and predictors response <- "Sale_Price" predictors <- setdiff(names(ames_data), response) ## ----------------------------------------------------------------------------- # Calculate correlations correlations <- sapply(predictors, function(predictor) { cor(ames_data[[predictor]], ames_data[[response]]) }) predictors <- names(correlations[abs(correlations) > 0.3]) ames_data <- ames_data[c(response, predictors)] print(names(ames_data)) ## ----------------------------------------------------------------------------- # Construct the formula for numeric predictors with function `f` applied predictors_formula_part <- paste0("f(", predictors, ")", collapse = " + ") # Combine both parts into the final formula ames_formula <- as.formula(paste(response, "~", predictors_formula_part)) print(ames_formula) ## ----message=FALSE------------------------------------------------------------ # Set the seed for reproducibility set.seed(123) # Determine the size of the dataset n <- nrow(ames_data) # Create a random sample of row indices for the training set trainIndex <- sample(1:n, size = floor(0.8 * n)) # Subset the data into training and test sets train <- ames_data[trainIndex, ] test <- ames_data[-trainIndex, ] Gmodboost <- NGeDSboost(formula = ames_formula, data = train, initial_learner = FALSE, phi_boost_exit = 0.999, phi = 0.95) ## ----------------------------------------------------------------------------- sqrt(mean((test[[response]] - predict(Gmodboost, newdata = test, n = 2))^2)) sqrt(mean((test[[response]] - predict(Gmodboost, newdata = test, n = 3))^2)) sqrt(mean((test[[response]] - predict(Gmodboost, newdata = test, n = 4))^2)) ## ----fig.width=10, fig.height=6----------------------------------------------- bl_imp <- bl_imp(Gmodboost) plot(bl_imp) ## ----fig.width=10, fig.height=6----------------------------------------------- # Linear FGB-GeDS Fit layout(matrix(c(1,2), nrow=1, byrow=TRUE)) plot(Gmodboost, n = 2) # Quadratic FGB-GeDS Fit plot(Gmodboost, n = 3) # Cubic FGB-GeDS Fit plot(Gmodboost, n = 4) ## ----------------------------------------------------------------------------- # Set the seed for reproducibility set.seed(123) # Determine the size of the dataset n <- nrow(ames_data) # Create a random sample of row indices for the training set trainIndex <- sample(1:n, size = floor(0.8 * n)) # Subset the data into training and test sets train <- ames_data[trainIndex, ] test <- ames_data[-trainIndex, ] Gmodboost <- NGeDSboost(formula = ames_formula, data = train, initial_learner = FALSE, phi = 0.95, shrinkage = 0.4) # RMSE sqrt(mean((test$Sale_Price - predict(Gmodboost, n = 2, newdata = test))^2)) ## ----------------------------------------------------------------------------- # Define order of the GeDS models ord = 3 # versicolor vs setosa: iris_subset <- subset(iris, Species %in% c("setosa", "versicolor")) iris_subset$Species <- factor(iris_subset$Species) mod <- glm(Species ~ Sepal.Length, family = binomial(link = "logit"), data = iris_subset) # GeDS Gmod <- suppressWarnings(GGeDS(Species ~ f(Sepal.Length), data = iris_subset, family = binomial(link = "logit"), phi = 0.8)) # GAM-GeDS Gmodgam <- suppressWarnings(NGeDSgam(Species ~ f(Sepal.Length), data = iris_subset, family = binomial(link = "logit"), phi = 0.8)) # FGB-GeDS Gmodboost <- suppressWarnings(NGeDSboost(Species ~ f(Sepal.Length), data = iris_subset, family = mboost::Binomial())) # Predicted multinomial probabilities pred_glm <- predict(mod, type = "response") pred_GeDS <- predict(Gmod, type = "response", n = ord) pred_GeDSgam <- predict(Gmodgam, type = "response", n = ord) pred_GeDSboost <- predict(Gmodboost, type = "response", n = ord) ## ----fig.width=10, fig.height=6----------------------------------------------- # Plot plot(iris_subset$Sepal.Length, as.numeric(iris_subset$Species) - 1, xlab = "Sepal Length", ylab = "Probability of Versicolor", main = "Classification: setosa vs. versicolor", pch = 16, col = "gray40") # GLM pred <- data.frame(Sepal.Length = as.vector(iris_subset$Sepal.Length), Predicted = pred_glm) pred <- pred[order(pred$Sepal.Length), ] lines(pred, col = "blue") # GeDS lines(pred$Sepal.Length, pred_GeDS, col = "red") # GAM-GeDS pred <- data.frame(Sepal.Length = Gmodboost$args$predictors$Sepal.Length, pred = pred_GeDSgam) pred <- pred[order(pred$Sepal.Length), ] lines(pred, col = "green") # FGB-GeDS pred <- data.frame(Sepal.Length = Gmodboost$args$predictors$Sepal.Length, ppred = pred_GeDSboost) pred <- pred[order(pred$Sepal.Length), ] lines(pred, col = "purple") # Legend legend("topleft", legend = c("GLM (logit)", "GeDS", "GAM-GeDS", "FGB-GeDS"), col = c("blue", "red", "green", "purple"), lwd = 2, bty = "n") ## ----------------------------------------------------------------------------- # Function to compute the confusion matrix and accuracy evaluate_model <- function(predicted, actual) { # Ensure both are factors with the same levels predicted <- factor(predicted, levels = levels(actual)) # Compute confusion matrix conf_matrix <- table(Predicted = predicted, Actual = actual) # Compute accuracy accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix) # Return results as a list return(list( Confusion_Matrix = conf_matrix, Accuracy = accuracy )) } ## ----------------------------------------------------------------------------- compute_multinomial_probs <- function(mod_1, mod_2, newdata, n) { # Compute linear predictors if (missing(n)) { eta1 <- predict(mod_1, newdata = newdata, type = "link") eta2 <- predict(mod_2, newdata = newdata, type = "link") } else { # For GeDS models the order needs to be specified eta1 <- predict(mod_1, newdata = newdata, type = "link", n = n) eta2 <- predict(mod_2, newdata = newdata, type = "link", n = n) } # Normalize exponentials for numerical stability (Log-Sum-Exp Trick) max_eta <- pmax(0, eta1, eta2) # Compute stabilized exponentials exp_setosa <- exp(0 - max_eta) # i.e., \approx 1 exp_versicolor <- exp(eta1 - max_eta) exp_virginica <- exp(eta2 - max_eta) # Calculate the denominator using the stabilized values denom <- exp_setosa + exp_versicolor + exp_virginica # Compute multinomial probabilities p_setosa <- exp_setosa / denom p_versicolor <- exp_versicolor / denom p_virginica <- exp_virginica / denom # Combine results in a data frame pred_probs <- data.frame( setosa = p_setosa, versicolor = p_versicolor, virginica = p_virginica ) return(pred_probs) } ## ----------------------------------------------------------------------------- multinomial_logit <- function(data, response, predictors, base_class, method = NGeDSboost, params = list()) { # Ensure response is a factor data[[response]] <- factor(data[[response]]) # Identify the other classes classes <- levels(data[[response]]) classes <- classes[classes != base_class] # Construct the formula dynamically with f() applied to each predictor predictor_formula <- paste(sprintf("f(%s)", predictors), collapse = " + ") formula_str <- paste(response, "~", predictor_formula) # Train binomial models models <- lapply(classes, function(class) { subset_data <- subset(data, data[[response]] %in% c(base_class, class)) subset_data[[response]] <- factor(subset_data[[response]]) suppressWarnings(do.call(method, c(list(as.formula(formula_str), data = subset_data), params))) }) # Compute multinomial probabilities pred_probs <- suppressWarnings(compute_multinomial_probs(models[[1]], models[[2]], newdata = data)) # Predict classes predicted_classes <- factor( apply(pred_probs, 1, function(x) colnames(pred_probs)[which.max(x)]), levels = colnames(pred_probs) ) # Evaluate model evaluation_results <- evaluate_model(predicted_classes, data[[response]]) return(list(probabilities = pred_probs, predicted_classes = predicted_classes, evaluation = evaluation_results)) } ## ----------------------------------------------------------------------------- # Example usage: # NGeDSboost results <- multinomial_logit(iris, response = "Species", predictors = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"), base_class = "setosa", method = NGeDSboost, params = list( family = mboost::Binomial(), phi = 0.99, q = 2 )) # Access results: print(head(results$probabilities)) print(head(results$predicted_classes)) print(results$evaluation)