diff --git a/DESCRIPTION b/DESCRIPTION
index ed36d62..b2b9c15 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: rwa
 Type: Package
 Title: Perform a Relative Weights Analysis
-Version: 0.1.1
+Version: 0.1.2.9000
 Authors@R: person(given = "Martin",
                   family = "Chan",
                   role = c("aut", "cre"),
diff --git a/NAMESPACE b/NAMESPACE
index 508b637..29c37a1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,8 @@ export("%>%")
 export(plot_rwa)
 export(remove_all_na_cols)
 export(rwa)
+export(rwa_logit)
+export(rwa_multiregress)
 import(dplyr)
 import(ggplot2)
 importFrom(boot,boot)
@@ -17,8 +19,14 @@ importFrom(dplyr,select)
 importFrom(dplyr,tibble)
 importFrom(magrittr,"%>%")
 importFrom(purrr,map_dfr)
+importFrom(stats,binomial)
+importFrom(stats,coef)
 importFrom(stats,cor)
+importFrom(stats,glm)
+importFrom(stats,lm)
+importFrom(stats,predict)
 importFrom(stats,rnorm)
+importFrom(stats,sd)
 importFrom(stats,var)
 importFrom(tidyr,drop_na)
 importFrom(utils,head)
diff --git a/NEWS.md b/NEWS.md
index 0873520..fd6a45f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,9 @@
+# rwa (development version)
+
+- Added `rwa_logit()` and `rwa_multiregress()` to support logistic regression and multiple regression. 
+- Added new vignette to cover the new regression methods. 
+- Improved test coverage and minor bugfixes.
+
 # rwa 0.1.1
 
 ## Improvements
diff --git a/R/globals.R b/R/globals.R
index 67bea59..ce9c9e2 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -9,6 +9,7 @@ utils::globalVariables(
     "Sign",
     "Rescaled.RelWeight",
     "Variables",
-    "Sign.Rescaled.RelWeight"
+    "Sign.Rescaled.RelWeight",
+    "variable_name"
   )
 )
diff --git a/R/plot_rwa.R b/R/plot_rwa.R
index de961b3..e630163 100644
--- a/R/plot_rwa.R
+++ b/R/plot_rwa.R
@@ -26,24 +26,48 @@ plot_rwa <- function(rwa){
     dplyr::mutate(Sign.Rescaled.RelWeight = ifelse(Sign == "-",
                                                    Rescaled.RelWeight * -1,
                                                    Rescaled.RelWeight))
+ 
+  # Get variable names from the Variables column
+  if (!"Variables" %in% names(result)) {
+    stop("Could not find 'Variables' column in result.")
+  }
+  result$variable_name <- result$Variables
 
+  # Calculate appropriate axis limits for both positive and negative values
+  max_abs_weight <- max(abs(result$Sign.Rescaled.RelWeight))
+  min_weight <- min(result$Sign.Rescaled.RelWeight)
   max_weight <- max(result$Sign.Rescaled.RelWeight)
+  
+  # Set symmetric or appropriate limits based on data
+  if (min_weight < 0 && max_weight > 0) {
+    # Mixed signs: use symmetric limits
+    y_limits <- c(-max_abs_weight * 1.15, max_abs_weight * 1.15)
+  } else if (min_weight < 0) {
+    # All negative: extend left for labels
+    y_limits <- c(min_weight * 1.15, max(0, max_weight * 1.1))
+  } else {
+    # All positive: extend right for labels
+    y_limits <- c(min(0, min_weight), max_weight * 1.15)
+  }
 
   result %>%
-    ggplot(aes(x = stats::reorder(Variables, Sign.Rescaled.RelWeight), 
+    ggplot(aes(x = stats::reorder(variable_name, Sign.Rescaled.RelWeight), 
                y = Sign.Rescaled.RelWeight)) +
     geom_col(fill = "#0066cc") +
-    geom_text(aes(label = round(Sign.Rescaled.RelWeight, 1)), hjust = -0.3) +
+    geom_text(aes(label = round(Sign.Rescaled.RelWeight, 1),
+                  hjust = ifelse(Sign.Rescaled.RelWeight >= 0, -0.3, 1.3))) +
     coord_flip() +
-    ylim(c(NA, max_weight * 1.1)) + # Automatic lower limit
+    ylim(y_limits) +
     labs(title = "Variable importance estimates",
          subtitle = "Using the Relative Weights Analysis method",
          x = "Predictor variables",
-         y = "Rescaled Relative Weights",
-         caption = paste0("Note: Rescaled Relative Weights sum to 100%. n = ",
+         y = "Rescaled Relative Weights (with sign)",
+         caption = paste0("Note: Absolute Rescaled Relative Weights sum to 100%. n = ",
                           rwa$n, ". ",
-                          "R-squared: ",
-                          round(rwa$rsquare, 2),
-                          ".")) +
+                          if (!is.null(rwa$rsquare)) {
+                            paste0("R-squared: ", round(rwa$rsquare, 2), ".")
+                          } else {
+                            ""
+                          })) +
     theme_classic()
 }
diff --git a/R/rwa.R b/R/rwa.R
index 875da4a..05c0f7a 100644
--- a/R/rwa.R
+++ b/R/rwa.R
@@ -1,51 +1,76 @@
 #' @title Create a Relative Weights Analysis (RWA)
 #'
-#' @description This function creates a Relative Weights Analysis (RWA) and 
-#'   returns a list of outputs. RWA provides a heuristic method for estimating 
-#'   the relative weight of predictor variables in multiple regression, which 
-#'   involves creating a multiple regression with on a set of transformed 
-#'   predictors which are orthogonal to each other but maximally related to the 
-#'   original set of predictors.
-#' `rwa()` is optimised for dplyr pipes and shows positive / negative signs for weights.
+#' @description This function creates a Relative Weights Analysis (RWA) and
+#'   returns a list of outputs. RWA provides a heuristic method for estimating
+#'   the relative weight of predictor variables in multiple regression, which
+#'   involves creating a multiple regression with on a set of transformed
+#'   predictors which are orthogonal to each other but maximally related to the
+#'   original set of predictors. `rwa()` is optimised for dplyr pipes and shows
+#'   positive / negative signs for weights.
 #'
 #' @details
-#' `rwa()` produces raw relative weight values (epsilons) as well as rescaled 
-#' weights (scaled as a percentage of predictable variance) for every predictor 
-#' in the model. Signs are added to the weights when the `applysigns` argument 
-#' is set to `TRUE`.
-#' See https://www.scotttonidandel.com/rwa-web for the 
+#' `rwa()` produces raw relative weight values (epsilons) as well as rescaled
+#' weights (scaled as a percentage of predictable variance) for every predictor
+#' in the model. Signs are added to the weights when the `applysigns` argument
+#' is set to `TRUE`. See <https://www.scotttonidandel.com/rwa-web> for the
 #' original implementation that inspired this package.
 #'
+#' This function is a wrapper around `rwa_multiregress()` and `rwa_logit()`,
+#' automatically selecting the appropriate method based on the outcome variable
+#' or the `method` argument.
+#'
 #' @param df Data frame or tibble to be passed through.
-#' @param outcome Outcome variable, to be specified as a string or bare input. 
+#' @param outcome Outcome variable, to be specified as a string or bare input.
 #'   Must be a numeric variable.
-#' @param predictors Predictor variable(s), to be specified as a vector of 
+#' @param predictors Predictor variable(s), to be specified as a vector of
 #'   string(s) or bare input(s). All variables must be numeric.
-#' @param applysigns Logical value specifying whether to show an estimate that applies the sign. Defaults to `FALSE`.
-#' @param sort Logical value specifying whether to sort results by rescaled relative weights in descending order. Defaults to `TRUE`.
-#' @param bootstrap Logical value specifying whether to calculate bootstrap confidence intervals. Defaults to `FALSE`.
-#' @param n_bootstrap Number of bootstrap samples to use when bootstrap = TRUE. Defaults to 1000.
+#' @param applysigns Logical value specifying whether to show an estimate that
+#'   applies the sign. Defaults to `FALSE`.
+#' @param method String to specify the method of regression to apply. Valid
+#'   values include:
+#'   - `"auto"`: automatically detect whether to use multiple regression or
+#'     logistic regression based on the outcome variable provided.
+#'   - `"multiple"`: use multiple regression.
+#'   - `"logistic"`: use logistic regression.
+#' @param sort Logical value specifying whether to sort results by rescaled
+#'   relative weights in descending order. Defaults to `TRUE`.
+#' @param bootstrap Logical value specifying whether to calculate bootstrap
+#'   confidence intervals. Defaults to `FALSE`. Currently only supported for
+#'   multiple regression.
+#' @param n_bootstrap Number of bootstrap samples to use when bootstrap = TRUE.
+#'   Defaults to 1000.
 #' @param conf_level Confidence level for bootstrap intervals. Defaults to 0.95.
 #' @param focal Focal variable for bootstrap comparisons (optional).
-#' @param comprehensive Whether to run comprehensive bootstrap analysis including random variable and focal comparisons.
-#' @param include_rescaled_ci Logical value specifying whether to include confidence intervals for rescaled weights. Defaults to `FALSE` due to compositional data constraints. Use with caution.
+#' @param comprehensive Whether to run comprehensive bootstrap analysis
+#'   including random variable and focal comparisons.
+#' @param include_rescaled_ci Logical value specifying whether to include
+#'   confidence intervals for rescaled weights. Defaults to `FALSE` due to
+#'   compositional data constraints. Use with caution.
 #'
 #' @return `rwa()` returns a list of outputs, as follows:
 #' - `predictors`: character vector of names of the predictor variables used.
-#' - `rsquare`: the rsquare value of the regression model.
-#' - `result`: the final output of the importance metrics (sorted by Rescaled.RelWeight in descending order by default).
+#' - `rsquare`: the rsquare value of the regression model (multiple regression only).
+#' - `result`: the final output of the importance metrics (sorted by
+#'   Rescaled.RelWeight in descending order by default).
 #'   - The `Rescaled.RelWeight` column sums up to 100.
-#'   - The `Sign` column indicates whether a predictor is positively or negatively correlated with the outcome.
+#'   - The `Sign` column indicates whether a predictor is positively or
+#'     negatively correlated with the outcome.
 #'   - When bootstrap = TRUE, includes confidence interval columns for raw weights.
-#'   - Rescaled weight CIs are available via include_rescaled_ci = TRUE but not recommended for inference.
+#'   - Rescaled weight CIs are available via include_rescaled_ci = TRUE but not
+#'     recommended for inference.
 #' - `n`: indicates the number of observations used in the analysis.
 #' - `bootstrap`: bootstrap results (only present when bootstrap = TRUE), containing:
 #'   - `ci_results`: confidence intervals for weights
 #'   - `boot_object`: raw bootstrap object for advanced analysis
 #'   - `n_bootstrap`: number of bootstrap samples used
-#' - `lambda`:
-#' - `RXX`: Correlation matrix of all the predictor variables against each other.
-#' - `RXY`: Correlation values of the predictor variables against the outcome variable.
+#' - `lambda`: lambda matrix from the RWA calculation.
+#' - `RXX`: Correlation matrix of all the predictor variables against each
+#'   other. Not available for logistic regression.
+#' - `RXY`: Correlation values of the predictor variables against the outcome
+#'   variable. Not available for logistic regression.
+#'
+#' @seealso [plot_rwa()] for plotting results, [rwa_multiregress()] and
+#'   [rwa_logit()] for the underlying implementations.
 #'
 #' @importFrom magrittr %>%
 #' @importFrom tidyr drop_na
@@ -54,180 +79,158 @@
 #' @examples
 #' library(ggplot2)
 #' # Basic RWA (results sorted by default)
-#' rwa(diamonds,"price",c("depth","carat"))
-#' 
+#' rwa(diamonds, "price", c("depth", "carat"))
+#'
 #' # RWA without sorting (preserves original predictor order)
-#' rwa(diamonds,"price",c("depth","carat"), sort = FALSE)
-#' 
+#' rwa(diamonds, "price", c("depth", "carat"), sort = FALSE)
+#'
+#' # Plot results using plot_rwa()
+#' diamonds |>
+#'   rwa("price", c("depth", "carat", "x", "y")) |>
+#'   plot_rwa()
+#'
 #' \donttest{
 #' # For faster examples, use a subset of data for bootstrap
 #' diamonds_small <- diamonds[sample(nrow(diamonds), 1000), ]
-#' 
+#'
 #' # RWA with bootstrap confidence intervals (raw weights only)
-#' rwa(diamonds_small,"price",c("depth","carat"), bootstrap = TRUE, n_bootstrap = 100)
-#' 
+#' rwa(diamonds_small, "price", c("depth", "carat"),
+#'     bootstrap = TRUE, n_bootstrap = 100)
+#'
 #' # Include rescaled weight CIs (use with caution for inference)
-#' rwa(diamonds_small,"price",c("depth","carat"), bootstrap = TRUE, 
-#'     include_rescaled_ci = TRUE, n_bootstrap = 100)
-#' 
+#' rwa(diamonds_small, "price", c("depth", "carat"),
+#'     bootstrap = TRUE, include_rescaled_ci = TRUE, n_bootstrap = 100)
+#'
 #' # Comprehensive bootstrap analysis with focal variable
-#' result <- rwa(diamonds_small,"price",c("depth","carat","table"), 
-#'               bootstrap = TRUE, comprehensive = TRUE, focal = "carat", 
+#' result <- rwa(diamonds_small, "price", c("depth", "carat", "table"),
+#'               bootstrap = TRUE, comprehensive = TRUE, focal = "carat",
 #'               n_bootstrap = 100)
 #' # View confidence intervals
 #' result$bootstrap$ci_results
 #' }
 #'
+#' # Based on logistic regression (auto-detected from binary outcome)
+#' diamonds$IsIdeal <- as.numeric(diamonds$cut == "Ideal")
+#' rwa(diamonds, "IsIdeal", c("depth", "carat"))
+#'
 #' @export
 rwa <- function(df,
                 outcome,
                 predictors,
                 applysigns = FALSE,
+                method = "auto",
                 sort = TRUE,
                 bootstrap = FALSE,
                 n_bootstrap = 1000,
                 conf_level = 0.95,
                 focal = NULL,
                 comprehensive = FALSE,
-                include_rescaled_ci = FALSE){
+                include_rescaled_ci = FALSE) {
 
 
   # ---- Input validation ----
-  
+
+  # Validate method
+
+  if (!method %in% c("auto", "multiple", "logistic")) {
+    stop("Invalid input for `method`. Must be one of: 'auto', 'multiple', 'logistic'.")
+  }
 
   # Validate conf_level
-  if (!is.numeric(conf_level) || length(conf_level) != 1 || 
+  if (!is.numeric(conf_level) || length(conf_level) != 1 ||
       conf_level <= 0 || conf_level >= 1) {
     stop("`conf_level` must be a single numeric value between 0 and 1 (exclusive).")
   }
-  
+
   # Validate n_bootstrap
-  if (!is.numeric(n_bootstrap) || length(n_bootstrap) != 1 || 
+  if (!is.numeric(n_bootstrap) || length(n_bootstrap) != 1 ||
       n_bootstrap < 1 || n_bootstrap != floor(n_bootstrap)) {
     stop("`n_bootstrap` must be a positive integer.")
   }
-  
+
   # Check that outcome and predictors exist in data
   if (!outcome %in% names(df)) {
     stop(sprintf("Outcome variable '%s' not found in data.", outcome))
   }
-  
+
   missing_predictors <- predictors[!predictors %in% names(df)]
   if (length(missing_predictors) > 0) {
-    stop(sprintf("Predictor variable(s) not found in data: %s", 
+    stop(sprintf("Predictor variable(s) not found in data: %s",
                  paste(missing_predictors, collapse = ", ")))
   }
-  
+
   # Validate that outcome is numeric
   if (!is.numeric(df[[outcome]])) {
     stop(sprintf("Outcome variable '%s' must be numeric.", outcome))
   }
-  
+
   # Validate that all predictors are numeric
   non_numeric_predictors <- predictors[!sapply(df[predictors], is.numeric)]
   if (length(non_numeric_predictors) > 0) {
-    stop(sprintf("All predictor variables must be numeric. Non-numeric: %s", 
+    stop(sprintf("All predictor variables must be numeric. Non-numeric: %s",
                  paste(non_numeric_predictors, collapse = ", ")))
   }
 
-  # Gets data frame in right order and form
-  thedata <-
-    df %>%
-    dplyr::select(dplyr::all_of(c(outcome, predictors))) %>%
-    tidyr::drop_na(dplyr::all_of(outcome))
+  # ---- Determine regression method ----
 
-  # Check for zero-variance outcome
-  outcome_var <- stats::var(thedata[[outcome]], na.rm = TRUE)
-  if (is.na(outcome_var) || outcome_var == 0) {
-    stop(sprintf("Outcome variable '%s' has zero variance.", outcome))
-  }
-  
-  # Check for zero-variance predictors
-  predictor_vars <- sapply(thedata[predictors], function(x) stats::var(x, na.rm = TRUE))
-  zero_var_predictors <- names(predictor_vars)[is.na(predictor_vars) | predictor_vars == 0]
-  if (length(zero_var_predictors) > 0) {
-    stop(sprintf("Predictor variable(s) with zero variance: %s", 
-                 paste(zero_var_predictors, collapse = ", ")))
-  }
-
-  cor_matrix <-
-    cor(thedata, use = "pairwise.complete.obs") %>%
-    as.data.frame(stringsAsFactors = FALSE, row.names = NULL) %>%
-    remove_all_na_cols() %>%
-    tidyr::drop_na()
-
-  matrix_data <-
-    cor_matrix %>%
-    as.matrix()
+  outcome_values <- unique(df[[outcome]])
+  outcome_var_unique <- dplyr::n_distinct(outcome_values[!is.na(outcome_values)])
 
-  RXX <- matrix_data[2:ncol(matrix_data), 2:ncol(matrix_data)] # Only take the correlations with the predictor variables
-  RXY <- matrix_data[2:ncol(matrix_data), 1] # Take the correlations of each of the predictors with the outcome variable
+  use_logistic <- (method == "logistic") ||
+                  (method == "auto" && outcome_var_unique == 2)
 
-  # Get all the 'genuine' predictor variables
-  Variables <-
-    cor_matrix %>%
-    names() %>%
-    .[.!=outcome]
-
-  # Check for singular/near-singular correlation matrix (perfect collinearity)
-  RXX_det <- det(as.matrix(RXX))
-  if (abs(RXX_det) < .Machine$double.eps * 100) {
-    stop("Predictor correlation matrix is singular or near-singular. This usually indicates perfect or near-perfect collinearity among predictors. Consider removing highly correlated predictors.")
-  }
-
-  RXX.eigen <- eigen(RXX) # Compute eigenvalues and eigenvectors of matrix
-  
-  # Check for negative eigenvalues (indicates numerical issues)
-  if (any(RXX.eigen$val < 0)) {
-    warning("Correlation matrix has negative eigenvalues, which may indicate numerical instability. Results should be interpreted with caution.")
+  if (method == "auto") {
+    if (use_logistic) {
+      message(
+        paste0("Parsing `", outcome, "` as a binary variable.\n"),
+        "Applying logistic regression to calculate relative weights..."
+      )
+    } else {
+      message(
+        paste0("Parsing `", outcome, "` as a non-binary variable.\n"),
+        "Applying multiple regression to calculate relative weights..."
+      )
+    }
   }
-  
-  D <- diag(RXX.eigen$val) # Run diag() on the values of eigen - construct diagonal matrix
-  delta <- sqrt(D) # Take square root of the created diagonal matrix
-
-  lambda <- RXX.eigen$vec %*% delta %*% t(RXX.eigen$vec) # Matrix multiplication
-  lambdasq <- lambda ^ 2 # Square the result
 
-  # To get partial effect of each independent variable on the dependent variable
-  # We multiply the inverse matrix (RXY) on the correlation matrix between dependent and independent variables
-  beta <- solve(lambda) %*% RXY # Solve numeric matrix containing coefficients of equation (Ax=B)
-  rsquare <- sum(beta ^ 2) # Output - R Square, sum of squared values
+  # ---- Handle bootstrap for logistic regression ----
 
-  RawWgt <- lambdasq %*% beta ^ 2 # Raw Relative Weight
-  import <- (RawWgt / rsquare) * 100 # Rescaled Relative Weight
+  if (bootstrap && use_logistic) {
+    warning("Bootstrap confidence intervals are not yet implemented for logistic regression. ",
+            "Proceeding without bootstrap.")
+    bootstrap <- FALSE
+  }
 
-  beta %>% # Get signs from coefficients
-    as.data.frame(stringsAsFactors = FALSE, row.names = NULL) %>%
-    dplyr::mutate_all(~(dplyr::case_when(.>0~"+",
-                                         .<0~"-",
-                                         .==0~"0",
-                           TRUE~NA_character_))) %>%
-    dplyr::rename(Sign="V1")-> sign
+  # ---- Call appropriate sub-function ----
 
-  result <- data.frame(Variables,
-                       Raw.RelWeight = RawWgt,
-                       Rescaled.RelWeight = import,
-                       Sign = sign) # Output - results
+  if (use_logistic) {
+    result_list <- rwa_logit(
+      df = df,
+      outcome = outcome,
+      predictors = predictors,
+      applysigns = applysigns
+    )
+  } else {
+    result_list <- rwa_multiregress(
+      df = df,
+      outcome = outcome,
+      predictors = predictors,
+      applysigns = applysigns
+    )
+  }
 
-  nrow(drop_na(thedata)) -> complete_cases
+  # ---- Apply sorting ----
 
-  if(applysigns == TRUE){
-    result %>%
-      dplyr::mutate(Sign.Rescaled.RelWeight = ifelse(Sign == "-",
-                                              Rescaled.RelWeight * -1,
-                                              Rescaled.RelWeight)) -> result
+  if (sort) {
+    result_list$result <- result_list$result %>%
+      dplyr::arrange(dplyr::desc(Rescaled.RelWeight))
   }
 
-  # Sort results by rescaled relative weights if requested
-  if(sort == TRUE){
-    result <- result %>%
-      dplyr::arrange(desc(Rescaled.RelWeight))
-  }
+  # ---- Run bootstrap analysis if requested (multiple regression only) ----
 
-  # Run bootstrap analysis if requested
-  if(bootstrap == TRUE) {
+  if (bootstrap) {
     message("Running bootstrap analysis with ", n_bootstrap, " samples...")
-    
+
     bootstrap_results <- run_rwa_bootstrap(
       data = df,
       outcome = outcome,
@@ -236,50 +239,41 @@ rwa <- function(df,
       conf_level = conf_level,
       focal = focal,
       comprehensive = comprehensive,
-      include_rescaled = include_rescaled_ci  # Only include if explicitly requested
+      include_rescaled = include_rescaled_ci
     )
-    
+
     # Add confidence intervals to result dataframe
-    if(!is.null(bootstrap_results$ci_results$raw_weights)) {
+    if (!is.null(bootstrap_results$ci_results$raw_weights)) {
       ci_data <- bootstrap_results$ci_results$raw_weights
-      
+
       # Add CI columns for raw weights
-      result$Raw.RelWeight.CI.Lower <- ci_data$ci_lower[match(result$Variables, ci_data$variable)]
-      result$Raw.RelWeight.CI.Upper <- ci_data$ci_upper[match(result$Variables, ci_data$variable)]
-      
+      result_list$result$Raw.RelWeight.CI.Lower <-
+        ci_data$ci_lower[match(result_list$result$Variables, ci_data$variable)]
+      result_list$result$Raw.RelWeight.CI.Upper <-
+        ci_data$ci_upper[match(result_list$result$Variables, ci_data$variable)]
+
       # Add significance indicator for raw weights (if CI doesn't include 0)
-      result$Raw.Significant <- !(result$Raw.RelWeight.CI.Lower <= 0 & result$Raw.RelWeight.CI.Upper >= 0)
+      result_list$result$Raw.Significant <-
+        !(result_list$result$Raw.RelWeight.CI.Lower <= 0 &
+          result_list$result$Raw.RelWeight.CI.Upper >= 0)
     }
-    
+
     # Add rescaled weight CIs only if explicitly requested and warn user
-    if(include_rescaled_ci && !is.null(bootstrap_results$ci_results$rescaled_weights)) {
-      warning("Rescaled weight confidence intervals should be interpreted with caution due to compositional data constraints. Use for descriptive purposes only, not formal statistical inference.")
-      
+    if (include_rescaled_ci && !is.null(bootstrap_results$ci_results$rescaled_weights)) {
+      warning("Rescaled weight confidence intervals should be interpreted with caution ",
+              "due to compositional data constraints. Use for descriptive purposes only, ",
+              "not formal statistical inference.")
+
       rescaled_ci <- bootstrap_results$ci_results$rescaled_weights
-      
-      result$Rescaled.RelWeight.CI.Lower <- rescaled_ci$ci_lower[match(result$Variables, rescaled_ci$variable)]
-      result$Rescaled.RelWeight.CI.Upper <- rescaled_ci$ci_upper[match(result$Variables, rescaled_ci$variable)]
-      
-      # Note: Not adding significance indicator for rescaled weights due to interpretation issues
+
+      result_list$result$Rescaled.RelWeight.CI.Lower <-
+        rescaled_ci$ci_lower[match(result_list$result$Variables, rescaled_ci$variable)]
+      result_list$result$Rescaled.RelWeight.CI.Upper <-
+        rescaled_ci$ci_upper[match(result_list$result$Variables, rescaled_ci$variable)]
     }
-    
-    return_list <- list("predictors" = Variables,
-                       "rsquare" = rsquare,
-                       "result" = result,
-                       "n" = complete_cases,
-                       "bootstrap" = bootstrap_results,
-                       "lambda" = lambda,
-                       "RXX" = RXX,
-                       "RXY" = RXY)
-  } else {
-    return_list <- list("predictors" = Variables,
-                       "rsquare" = rsquare,
-                       "result" = result,
-                       "n" = complete_cases,
-                       "lambda" = lambda,
-                       "RXX" = RXX,
-                       "RXY" = RXY)
+
+    result_list$bootstrap <- bootstrap_results
   }
 
-  return_list
-}
\ No newline at end of file
+  result_list
+}
diff --git a/R/rwa_logit.R b/R/rwa_logit.R
new file mode 100644
index 0000000..aeb02a2
--- /dev/null
+++ b/R/rwa_logit.R
@@ -0,0 +1,165 @@
+#' @title Create a Relative Weights Analysis with logistic regression
+#'
+#' @description This function performs Relative Weights Analysis (RWA) for binary
+#'   outcome variables using logistic regression. RWA provides a method for
+#'   estimating the relative importance of predictor variables by transforming
+#'   them into orthogonal variables while preserving their relationship to the
+#'   outcome. This implementation follows Johnson (2000) for logistic regression.
+#'
+#' @inheritParams rwa
+#'
+#' @return `rwa_logit()` returns a list of outputs, as follows:
+#' - `predictors`: character vector of names of the predictor variables used.
+#' - `rsquare`: the pseudo R-squared value (sum of epsilon weights) for the logistic regression model.
+#' - `result`: the final output of the importance metrics.
+#'   - The `Rescaled.RelWeight` column sums up to 100.
+#'   - The `Sign` column indicates whether a predictor is positively or negatively associated with the outcome.
+#' - `n`: indicates the number of observations used in the analysis.
+#' - `lambda`: the Lambda transformation matrix from the analysis.
+#'
+#' @examples
+#' # Create a binary outcome variable
+#' mtcars_binary <- mtcars
+#' mtcars_binary$high_mpg <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)
+#'
+#' # Basic logistic RWA
+#' result <- rwa_logit(
+#'   df = mtcars_binary,
+#'   outcome = "high_mpg",
+#'   predictors = c("cyl", "disp", "hp", "wt")
+#' )
+#'
+#' # View the relative importance results
+#' result$result
+#'
+#' # With sign information
+#' result_signed <- rwa_logit(
+#'   df = mtcars_binary,
+#'   outcome = "high_mpg",
+#'   predictors = c("cyl", "disp", "hp", "wt"),
+#'   applysigns = TRUE
+#' )
+#' result_signed$result
+#'
+#' @importFrom magrittr %>%
+#' @importFrom stats glm binomial coef predict sd lm
+#' @export
+rwa_logit <- function(df,
+                      outcome,
+                      predictors,
+                      applysigns = FALSE){
+
+  # Gets data frame in right order and form
+  thedata <-
+    df %>%
+    dplyr::select(all_of(c(outcome, predictors))) %>%
+    tidyr::drop_na(all_of(outcome))
+
+  # Get variable names for output
+  Variables <-
+    thedata %>%
+    dplyr::select(all_of(predictors)) %>%
+    names()
+
+  # Select outcome variable
+  Y <-
+    thedata %>%
+    pull(all_of(outcome))
+
+  # Scaled predictors
+  X <-
+    thedata %>%
+    dplyr::select(all_of(predictors)) %>%
+    scale()
+
+  X.svd <- svd(X) # Single-value decomposition
+  Q <- X.svd$v
+  P <- X.svd$u
+  Z <- P %*% t(Q)
+
+  Z.stand <- scale(Z)
+
+  # Obtaining Lambda from equation 7 from Johnson (2000) pg 8
+  Lambda <-
+    solve(
+      t(Z.stand) %*% Z.stand
+      ) %*%
+    t(Z.stand) %*%
+    X
+
+  logrfit <-
+    glm(Y ~ Z.stand,
+        family = binomial)
+
+  unstCoefs <- coef(logrfit)
+
+  b <- unstCoefs[2:length(unstCoefs)]
+
+  LpredY <-
+    predict(
+      logrfit,
+      newdata = thedata,
+      type="response")
+
+  # Clamp predictions to avoid Inf/-Inf in logit transformation
+  # This can occur with perfect separation in logistic regression
+  LpredY <- pmax(pmin(LpredY, 1 - 1e-10), 1e-10)
+
+  # Creating logit-Y-hat
+  lYhat <- log(LpredY/(1-LpredY))
+
+  # Getting st dev of logit-Y-hat
+  stdlYhat <- stats::sd(lYhat)
+
+  # Check for zero standard deviation (can occur with perfect separation)
+  if (stdlYhat == 0 || is.na(stdlYhat)) {
+    warning("Perfect or near-perfect separation detected. Results may be unreliable.")
+    stdlYhat <- 1e-10  # Avoid division by zero
+  }
+
+  # Getting R-sq
+  getting.Rsq <- lm(LpredY ~ Y)
+
+  # Computing standardized logistic regression coefficients
+  Rsq <- summary(getting.Rsq)$r.squared
+  beta <- b*((sqrt(Rsq))/stdlYhat)
+
+  epsilon <- Lambda^2 %*% beta^2
+
+  R.sq <- sum(epsilon)
+  PropWeights <- (epsilon/R.sq) * 100  # Convert to percentage (0-100) for consistency with rwa_multiregress
+
+  # Get signs from coefficients
+  sign <-
+    data.frame(V1 = beta) %>%
+    dplyr::mutate_all(~(dplyr::case_when(.>0~"+",
+                                         .<0~"-",
+                                         .==0~"0",
+                                         TRUE~NA_character_))) %>%
+    dplyr::rename(Sign = "V1")
+
+  ## Result
+  result <-
+    data.frame(Variables,
+               Raw.RelWeight = epsilon,
+               Rescaled.RelWeight = PropWeights) %>%
+    dplyr::mutate(Sign = sign)
+
+  complete_cases <- nrow(tidyr::drop_na(thedata))
+
+  if(applysigns == TRUE){
+    result <-
+      result %>%
+      dplyr::mutate(Sign.Rescaled.RelWeight = ifelse(Sign == "-",
+                                              Rescaled.RelWeight * -1,
+                                              Rescaled.RelWeight))
+  }
+
+  lambda <- Lambda
+
+  list("predictors" = predictors,
+       "rsquare" = R.sq,
+       "result" = result,
+       "n" = complete_cases,
+       "lambda" = lambda)
+}
diff --git a/R/rwa_multiregress.R b/R/rwa_multiregress.R
new file mode 100644
index 0000000..4f2f108
--- /dev/null
+++ b/R/rwa_multiregress.R
@@ -0,0 +1,131 @@
+#' @title Create a Relative Weights Analysis (RWA)
+#'
+#' @description This function creates a Relative Weights Analysis (RWA) and returns a list of outputs.
+#' RWA provides a heuristic method for estimating the relative weight of predictor variables in multiple regression, which involves
+#' creating a multiple regression with on a set of transformed predictors which are orthogonal to each other but
+#' maximally related to the original set of predictors.
+#' `rwa_multiregress()` is optimised for dplyr pipes and shows positive / negative signs for weights.
+#'
+#' @details
+#' `rwa_multiregress()` produces raw relative weight values (epsilons) as well as rescaled weights (scaled as a percentage of predictable variance)
+#' for every predictor in the model.
+#' Signs are added to the weights when the `applysigns` argument is set to `TRUE`.
+#' See https://relativeimportance.davidson.edu/multipleregression.html for the original implementation that inspired this package.
+#'
+#' @param df Data frame or tibble to be passed through.
+#' @param outcome Outcome variable, to be specified as a string or bare input. Must be a numeric variable.
+#' @param predictors Predictor variable(s), to be specified as a vector of string(s) or bare input(s). All variables must be numeric.
+#' @param applysigns Logical value specifying whether to show an estimate that applies the sign. Defaults to `FALSE`.
+#'
+#' @return `rwa_multiregress()` returns a list of outputs, as follows:
+#' - `predictors`: character vector of names of the predictor variables used.
+#' - `rsquare`: the rsquare value of the regression model.
+#' - `result`: the final output of the importance metrics.
+#'   - The `Rescaled.RelWeight` column sums up to 100.
+#'   - The `Sign` column indicates whether a predictor is positively or negatively correlated with the outcome.
+#' - `n`: indicates the number of observations used in the analysis.
+#' - `lambda`: the transformation matrix that maps the original correlated predictors to orthogonal variables while preserving their relationship to the outcome. Used internally to compute relative weights.
+#' - `RXX`: Correlation matrix of all the predictor variables against each other.
+#' - `RXY`: Correlation values of the predictor variables against the outcome variable.
+#'
+#' @importFrom magrittr %>%
+#' @importFrom tidyr drop_na
+#' @importFrom stats cor
+#' @import dplyr
+#' @examples
+#' # Basic multiple regression RWA
+#' result <- rwa_multiregress(
+#'   df = mtcars,
+#'   outcome = "mpg",
+#'   predictors = c("cyl", "disp", "hp", "wt")
+#' )
+#'
+#' # View the relative importance results
+#' result$result
+#'
+#' # With sign information
+#' result_signed <- rwa_multiregress(
+#'   df = mtcars,
+#'   outcome = "mpg",
+#'   predictors = c("cyl", "disp", "hp", "wt"),
+#'   applysigns = TRUE
+#' )
+#' result_signed$result
+#'
+#' @export
+rwa_multiregress <- function(df,
+                             outcome,
+                             predictors,
+                             applysigns = FALSE){
+
+  # Gets data frame in right order and form
+  thedata <-
+    df %>%
+    dplyr::select(all_of(c(outcome, predictors))) %>%
+    tidyr::drop_na(all_of(outcome))
+
+  cor_matrix <-
+    cor(thedata, use = "pairwise.complete.obs") %>%
+    as.data.frame(stringsAsFactors = FALSE, row.names = NULL) %>%
+    remove_all_na_cols() %>%
+    tidyr::drop_na()
+
+  matrix_data <-
+    cor_matrix %>%
+    as.matrix()
+
+  RXX <- matrix_data[2:ncol(matrix_data), 2:ncol(matrix_data)] # Only take the correlations with the predictor variables
+  RXY <- matrix_data[2:ncol(matrix_data), 1] # Take the correlations of each of the predictors with the outcome variable
+
+  # Get all the 'genuine' predictor variables
+  Variables <-
+    cor_matrix %>%
+    names() %>%
+    .[.!=outcome]
+
+  RXX.eigen <- eigen(RXX) # Compute eigenvalues and eigenvectors of matrix
+  D <- diag(RXX.eigen$val) # Run diag() on the values of eigen - construct diagonal matrix
+  delta <- sqrt(D) # Take square root of the created diagonal matrix
+
+  lambda <- RXX.eigen$vec %*% delta %*% t(RXX.eigen$vec) # Matrix multiplication
+  lambdasq <- lambda ^ 2 # Square the result
+
+  # To get partial effect of each independent variable on the dependent variable
+  # We multiply the inverse matrix (RXY) on the correlation matrix between dependent and independent variables
+  beta <- solve(lambda) %*% RXY # Solve numeric matrix containing coefficients of equation (Ax=B)
+  rsquare <- sum(beta ^ 2) # Output - R Square, sum of squared values
+
+  RawWgt <- lambdasq %*% beta ^ 2 # Raw Relative Weight
+  import <- (RawWgt / rsquare) * 100 # Rescaled Relative Weight
+
+  sign <- beta %>% # Get signs from coefficients
+    as.data.frame(stringsAsFactors = FALSE, row.names = NULL) %>%
+    dplyr::mutate_all(~(dplyr::case_when(.>0~"+",
+                                         .<0~"-",
+                                         .==0~"0",
+                           TRUE~NA_character_))) %>%
+    dplyr::rename(Sign="V1")
+
+  result <- data.frame(Variables,
+                       Raw.RelWeight = RawWgt,
+                       Rescaled.RelWeight = import,
+                       Sign = sign) # Output - results
+
+  complete_cases <- nrow(tidyr::drop_na(thedata))
+
+  if(applysigns == TRUE){
+    result <-
+      result %>%
+      dplyr::mutate(Sign.Rescaled.RelWeight = ifelse(Sign == "-",
+                                              Rescaled.RelWeight * -1,
+                                              Rescaled.RelWeight))
+  }
+
+  list("predictors" = Variables,
+       "rsquare" = rsquare,
+       "result" = result,
+       "n" = complete_cases,
+       "lambda" = lambda,
+       "RXX" = RXX,
+       "RXY" = RXY)
+}
diff --git a/README.md b/README.md
index 61a1079..2d310a3 100644
--- a/README.md
+++ b/README.md
@@ -152,5 +152,3 @@ Tonidandel, S., LeBreton, J. M., & Johnson, J. W. (2009). Determining the statis
 
 Wang, X., Duverger, P., Bansal, H. S. (2013). Bayesian Inference of Predictors Relative Importance in Linear Regression Model using Dominance Hierarchies. *International Journal of Pure and Applied Mathematics*, Vol. 88, No. 3, 321-339.
 
-Also see [Kovalyshyn](https://medium.com/analytics-vidhya/johnsons-relative-weights-analysis-implementation-with-javascript-d85393c0bbb4) for a similar implementation but in Javascript.
-
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 5ac9da4..e04d585 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -15,6 +15,8 @@ reference:
     desc: "Core functions for relative weights analysis"
     contents:
     - rwa
+    - rwa_multiregress
+    - rwa_logit
     - plot_rwa
   - title: "Utility Functions"
     desc: "Helper functions and utilities"
@@ -28,5 +30,6 @@ articles:
   
 - title: "Advanced Methods"
   contents:
+  - regression-methods
   - bootstrap-confidence-intervals
   - evaluating-rwa-method-reference
diff --git a/man/rwa.Rd b/man/rwa.Rd
index 8acccb5..a01aa8c 100644
--- a/man/rwa.Rd
+++ b/man/rwa.Rd
@@ -9,6 +9,7 @@ rwa(
   outcome,
   predictors,
   applysigns = FALSE,
+  method = "auto",
   sort = TRUE,
   bootstrap = FALSE,
   n_bootstrap = 1000,
@@ -27,33 +28,53 @@ Must be a numeric variable.}
 \item{predictors}{Predictor variable(s), to be specified as a vector of
 string(s) or bare input(s). All variables must be numeric.}
 
-\item{applysigns}{Logical value specifying whether to show an estimate that applies the sign. Defaults to \code{FALSE}.}
+\item{applysigns}{Logical value specifying whether to show an estimate that
+applies the sign. Defaults to \code{FALSE}.}
 
-\item{sort}{Logical value specifying whether to sort results by rescaled relative weights in descending order. Defaults to \code{TRUE}.}
+\item{method}{String to specify the method of regression to apply. Valid
+values include:
+\itemize{
+\item \code{"auto"}: automatically detect whether to use multiple regression or
+logistic regression based on the outcome variable provided.
+\item \code{"multiple"}: use multiple regression.
+\item \code{"logistic"}: use logistic regression.
+}}
+
+\item{sort}{Logical value specifying whether to sort results by rescaled
+relative weights in descending order. Defaults to \code{TRUE}.}
 
-\item{bootstrap}{Logical value specifying whether to calculate bootstrap confidence intervals. Defaults to \code{FALSE}.}
+\item{bootstrap}{Logical value specifying whether to calculate bootstrap
+confidence intervals. Defaults to \code{FALSE}. Currently only supported for
+multiple regression.}
 
-\item{n_bootstrap}{Number of bootstrap samples to use when bootstrap = TRUE. Defaults to 1000.}
+\item{n_bootstrap}{Number of bootstrap samples to use when bootstrap = TRUE.
+Defaults to 1000.}
 
 \item{conf_level}{Confidence level for bootstrap intervals. Defaults to 0.95.}
 
 \item{focal}{Focal variable for bootstrap comparisons (optional).}
 
-\item{comprehensive}{Whether to run comprehensive bootstrap analysis including random variable and focal comparisons.}
+\item{comprehensive}{Whether to run comprehensive bootstrap analysis
+including random variable and focal comparisons.}
 
-\item{include_rescaled_ci}{Logical value specifying whether to include confidence intervals for rescaled weights. Defaults to \code{FALSE} due to compositional data constraints. Use with caution.}
+\item{include_rescaled_ci}{Logical value specifying whether to include
+confidence intervals for rescaled weights. Defaults to \code{FALSE} due to
+compositional data constraints. Use with caution.}
 }
 \value{
 \code{rwa()} returns a list of outputs, as follows:
 \itemize{
 \item \code{predictors}: character vector of names of the predictor variables used.
-\item \code{rsquare}: the rsquare value of the regression model.
-\item \code{result}: the final output of the importance metrics (sorted by Rescaled.RelWeight in descending order by default).
+\item \code{rsquare}: the rsquare value of the regression model (multiple regression only).
+\item \code{result}: the final output of the importance metrics (sorted by
+Rescaled.RelWeight in descending order by default).
 \itemize{
 \item The \code{Rescaled.RelWeight} column sums up to 100.
-\item The \code{Sign} column indicates whether a predictor is positively or negatively correlated with the outcome.
+\item The \code{Sign} column indicates whether a predictor is positively or
+negatively correlated with the outcome.
 \item When bootstrap = TRUE, includes confidence interval columns for raw weights.
-\item Rescaled weight CIs are available via include_rescaled_ci = TRUE but not recommended for inference.
+\item Rescaled weight CIs are available via include_rescaled_ci = TRUE but not
+recommended for inference.
 }
 \item \code{n}: indicates the number of observations used in the analysis.
 \item \code{bootstrap}: bootstrap results (only present when bootstrap = TRUE), containing:
@@ -62,9 +83,11 @@ string(s) or bare input(s). All variables must be numeric.}
 \item \code{boot_object}: raw bootstrap object for advanced analysis
 \item \code{n_bootstrap}: number of bootstrap samples used
 }
-\item \code{lambda}:
-\item \code{RXX}: Correlation matrix of all the predictor variables against each other.
-\item \code{RXY}: Correlation values of the predictor variables against the outcome variable.
+\item \code{lambda}: lambda matrix from the RWA calculation.
+\item \code{RXX}: Correlation matrix of all the predictor variables against each
+other. Not available for logistic regression.
+\item \code{RXY}: Correlation values of the predictor variables against the outcome
+variable. Not available for logistic regression.
 }
 }
 \description{
@@ -73,42 +96,59 @@ returns a list of outputs. RWA provides a heuristic method for estimating
 the relative weight of predictor variables in multiple regression, which
 involves creating a multiple regression with on a set of transformed
 predictors which are orthogonal to each other but maximally related to the
-original set of predictors.
-\code{rwa()} is optimised for dplyr pipes and shows positive / negative signs for weights.
+original set of predictors. \code{rwa()} is optimised for dplyr pipes and shows
+positive / negative signs for weights.
 }
 \details{
 \code{rwa()} produces raw relative weight values (epsilons) as well as rescaled
 weights (scaled as a percentage of predictable variance) for every predictor
 in the model. Signs are added to the weights when the \code{applysigns} argument
-is set to \code{TRUE}.
-See https://www.scotttonidandel.com/rwa-web for the
+is set to \code{TRUE}. See \url{https://www.scotttonidandel.com/rwa-web} for the
 original implementation that inspired this package.
+
+This function is a wrapper around \code{rwa_multiregress()} and \code{rwa_logit()},
+automatically selecting the appropriate method based on the outcome variable
+or the \code{method} argument.
 }
 \examples{
 library(ggplot2)
 # Basic RWA (results sorted by default)
-rwa(diamonds,"price",c("depth","carat"))
+rwa(diamonds, "price", c("depth", "carat"))
 
 # RWA without sorting (preserves original predictor order)
-rwa(diamonds,"price",c("depth","carat"), sort = FALSE)
+rwa(diamonds, "price", c("depth", "carat"), sort = FALSE)
+
+# Plot results using plot_rwa()
+diamonds |>
+  rwa("price", c("depth", "carat", "x", "y")) |>
+  plot_rwa()
 
 \donttest{
 # For faster examples, use a subset of data for bootstrap
 diamonds_small <- diamonds[sample(nrow(diamonds), 1000), ]
 
 # RWA with bootstrap confidence intervals (raw weights only)
-rwa(diamonds_small,"price",c("depth","carat"), bootstrap = TRUE, n_bootstrap = 100)
+rwa(diamonds_small, "price", c("depth", "carat"),
+    bootstrap = TRUE, n_bootstrap = 100)
 
 # Include rescaled weight CIs (use with caution for inference)
-rwa(diamonds_small,"price",c("depth","carat"), bootstrap = TRUE, 
-    include_rescaled_ci = TRUE, n_bootstrap = 100)
+rwa(diamonds_small, "price", c("depth", "carat"),
+    bootstrap = TRUE, include_rescaled_ci = TRUE, n_bootstrap = 100)
 
 # Comprehensive bootstrap analysis with focal variable
-result <- rwa(diamonds_small,"price",c("depth","carat","table"), 
-              bootstrap = TRUE, comprehensive = TRUE, focal = "carat", 
+result <- rwa(diamonds_small, "price", c("depth", "carat", "table"),
+              bootstrap = TRUE, comprehensive = TRUE, focal = "carat",
               n_bootstrap = 100)
 # View confidence intervals
 result$bootstrap$ci_results
 }
 
+# Based on logistic regression (auto-detected from binary outcome)
+diamonds$IsIdeal <- as.numeric(diamonds$cut == "Ideal")
+rwa(diamonds, "IsIdeal", c("depth", "carat"))
+
+}
+\seealso{
+\code{\link[=plot_rwa]{plot_rwa()}} for plotting results, \code{\link[=rwa_multiregress]{rwa_multiregress()}} and
+\code{\link[=rwa_logit]{rwa_logit()}} for the underlying implementations.
 }
diff --git a/man/rwa_logit.Rd b/man/rwa_logit.Rd
new file mode 100644
index 0000000..6e1b411
--- /dev/null
+++ b/man/rwa_logit.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rwa_logit.R
+\name{rwa_logit}
+\alias{rwa_logit}
+\title{Create a Relative Weights Analysis with logistic regression}
+\usage{
+rwa_logit(df, outcome, predictors, applysigns = FALSE)
+}
+\arguments{
+\item{df}{Data frame or tibble to be passed through.}
+
+\item{outcome}{Outcome variable, to be specified as a string or bare input.
+Must be a numeric variable.}
+
+\item{predictors}{Predictor variable(s), to be specified as a vector of
+string(s) or bare input(s). All variables must be numeric.}
+
+\item{applysigns}{Logical value specifying whether to show an estimate that
+applies the sign. Defaults to \code{FALSE}.}
+}
+\value{
+\code{rwa_logit()} returns a list of outputs, as follows:
+\itemize{
+\item \code{predictors}: character vector of names of the predictor variables used.
+\item \code{rsquare}: the pseudo R-squared value (sum of epsilon weights) for the logistic regression model.
+\item \code{result}: the final output of the importance metrics.
+\itemize{
+\item The \code{Rescaled.RelWeight} column sums up to 100.
+\item The \code{Sign} column indicates whether a predictor is positively or negatively associated with the outcome.
+}
+\item \code{n}: indicates the number of observations used in the analysis.
+\item \code{lambda}: the Lambda transformation matrix from the analysis.
+}
+}
+\description{
+This function performs Relative Weights Analysis (RWA) for binary
+outcome variables using logistic regression. RWA provides a method for
+estimating the relative importance of predictor variables by transforming
+them into orthogonal variables while preserving their relationship to the
+outcome. This implementation follows Johnson (2000) for logistic regression.
+}
+\examples{
+# Create a binary outcome variable
+mtcars_binary <- mtcars
+mtcars_binary$high_mpg <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)
+
+# Basic logistic RWA
+result <- rwa_logit(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+
+# View the relative importance results
+result$result
+
+# With sign information
+result_signed <- rwa_logit(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  applysigns = TRUE
+)
+result_signed$result
+
+}
diff --git a/man/rwa_multiregress.Rd b/man/rwa_multiregress.Rd
new file mode 100644
index 0000000..da3353a
--- /dev/null
+++ b/man/rwa_multiregress.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rwa_multiregress.R
+\name{rwa_multiregress}
+\alias{rwa_multiregress}
+\title{Create a Relative Weights Analysis (RWA)}
+\usage{
+rwa_multiregress(df, outcome, predictors, applysigns = FALSE)
+}
+\arguments{
+\item{df}{Data frame or tibble to be passed through.}
+
+\item{outcome}{Outcome variable, to be specified as a string or bare input. Must be a numeric variable.}
+
+\item{predictors}{Predictor variable(s), to be specified as a vector of string(s) or bare input(s). All variables must be numeric.}
+
+\item{applysigns}{Logical value specifying whether to show an estimate that applies the sign. Defaults to \code{FALSE}.}
+}
+\value{
+\code{rwa_multiregress()} returns a list of outputs, as follows:
+\itemize{
+\item \code{predictors}: character vector of names of the predictor variables used.
+\item \code{rsquare}: the rsquare value of the regression model.
+\item \code{result}: the final output of the importance metrics.
+\itemize{
+\item The \code{Rescaled.RelWeight} column sums up to 100.
+\item The \code{Sign} column indicates whether a predictor is positively or negatively correlated with the outcome.
+}
+\item \code{n}: indicates the number of observations used in the analysis.
+\item \code{lambda}: the transformation matrix that maps the original correlated predictors to orthogonal variables while preserving their relationship to the outcome. Used internally to compute relative weights.
+\item \code{RXX}: Correlation matrix of all the predictor variables against each other.
+\item \code{RXY}: Correlation values of the predictor variables against the outcome variable.
+}
+}
+\description{
+This function creates a Relative Weights Analysis (RWA) and returns a list of outputs.
+RWA provides a heuristic method for estimating the relative weight of predictor variables in multiple regression, which involves
+creating a multiple regression with on a set of transformed predictors which are orthogonal to each other but
+maximally related to the original set of predictors.
+\code{rwa_multiregress()} is optimised for dplyr pipes and shows positive / negative signs for weights.
+}
+\details{
+\code{rwa_multiregress()} produces raw relative weight values (epsilons) as well as rescaled weights (scaled as a percentage of predictable variance)
+for every predictor in the model.
+Signs are added to the weights when the \code{applysigns} argument is set to \code{TRUE}.
+See https://relativeimportance.davidson.edu/multipleregression.html for the original implementation that inspired this package.
+}
+\examples{
+# Basic multiple regression RWA
+result <- rwa_multiregress(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+
+# View the relative importance results
+result$result
+
+# With sign information
+result_signed <- rwa_multiregress(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  applysigns = TRUE
+)
+result_signed$result
+
+}
diff --git a/tests/testthat/test-plot_rwa.R b/tests/testthat/test-plot_rwa.R
index 72231c7..76b7955 100644
--- a/tests/testthat/test-plot_rwa.R
+++ b/tests/testthat/test-plot_rwa.R
@@ -55,14 +55,14 @@ test_that("plot_rwa() plot structure is correct", {
   # Check plot data
   plot_data <- p$data
   expect_s3_class(plot_data, "data.frame")
-  expect_true("Variables" %in% names(plot_data))
+  expect_true("variable_name" %in% names(plot_data))
   expect_true("Sign.Rescaled.RelWeight" %in% names(plot_data))
   
   # Check that plot has correct number of bars
   expect_equal(nrow(plot_data), length(rwa_result$predictors))
   
   # Check mapping using rlang::as_label instead of as.character
-  expect_equal(rlang::as_label(p$mapping$x), "stats::reorder(Variables, Sign.Rescaled.RelWeight)")
+  expect_equal(rlang::as_label(p$mapping$x), "stats::reorder(variable_name, Sign.Rescaled.RelWeight)")
   expect_equal(rlang::as_label(p$mapping$y), "Sign.Rescaled.RelWeight")
   
   # Check coordinate system (should be flipped)
@@ -102,10 +102,10 @@ test_that("plot_rwa() labels and aesthetics are correct", {
   expect_equal(p$labels$title, "Variable importance estimates")
   expect_equal(p$labels$subtitle, "Using the Relative Weights Analysis method")
   expect_equal(p$labels$x, "Predictor variables")
-  expect_equal(p$labels$y, "Rescaled Relative Weights")
+  expect_equal(p$labels$y, "Rescaled Relative Weights (with sign)")
   
   # Check that caption contains expected elements
-  expect_true(grepl("Rescaled Relative Weights sum to 100%", p$labels$caption))
+  expect_true(grepl("Absolute Rescaled Relative Weights sum to 100%", p$labels$caption))
   expect_true(grepl("n =", p$labels$caption))
   expect_true(grepl("R-squared:", p$labels$caption))
   
@@ -146,15 +146,20 @@ test_that("plot_rwa() y-axis limits work correctly", {
   plot_data <- p$data
   
   # Check that y-limits accommodate the data plus some buffer
+  min_weight <- min(plot_data$Sign.Rescaled.RelWeight)
   max_weight <- max(plot_data$Sign.Rescaled.RelWeight)
   
-  # The ylim should be from NA (automatic) to max_weight * 1.1
-  # But if max_weight is negative (all negative relationships), 
-  # we should test against the absolute maximum for proper scaling
-  expected_upper_limit <- max_weight * 1.1
+  # The ylim should properly accommodate all values with buffer for labels
   y_limits <- layer_scales(p)$y$limits
   if (!is.null(y_limits)) {
-    expect_equal(y_limits[2], expected_upper_limit, tolerance = 1e-10)
+    # For all-negative case, lower limit should extend past minimum
+    if (min_weight < 0 && max_weight <= 0) {
+      expect_true(y_limits[1] < min_weight)
+    }
+    # For all-positive case, upper limit should extend past maximum
+    if (min_weight >= 0 && max_weight > 0) {
+      expect_true(y_limits[2] > max_weight)
+    }
   }
 })
 
@@ -173,3 +178,71 @@ test_that("plot_rwa() text labels are present", {
   text_layer <- p$layers[geom_classes == "GeomText"][[1]]
   expect_true("label" %in% names(text_layer$mapping))
 })
+
+# Tests for logistic regression plotting
+test_that("plot_rwa() works with logistic regression results", {
+  # Create binary outcome
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)
+  
+  # Run logistic RWA
+  rwa_logit_result <- rwa_logit(
+    df = mtcars_binary,
+    outcome = "high_mpg",
+    predictors = c("cyl", "hp", "wt")
+  )
+  
+  # plot_rwa() should work without error
+  expect_no_error(p <- plot_rwa(rwa_logit_result))
+  expect_s3_class(p, "ggplot")
+  
+  # Check that all predictors are represented
+  expect_equal(nrow(p$data), 3)
+})
+
+test_that("plot_rwa() caption shows R-squared for logistic regression", {
+  # Create binary outcome
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)
+  
+  # Run logistic RWA
+  rwa_logit_result <- rwa_logit(
+    df = mtcars_binary,
+    outcome = "high_mpg",
+    predictors = c("cyl", "hp", "wt")
+  )
+  
+  p <- plot_rwa(rwa_logit_result)
+  
+  # Caption should contain R-squared (logistic regression now returns rsquare too)
+  expect_true(grepl("R-squared:", p$labels$caption))
+})
+
+test_that("plot_rwa() caption shows R-squared for multiple regression", {
+  rwa_result <- mtcars %>%
+    rwa(outcome = "mpg", predictors = c("cyl", "hp", "wt"))
+  
+  p <- plot_rwa(rwa_result)
+  
+  # Caption should contain R-squared, not Lambda
+  expect_true(grepl("R-squared:", p$labels$caption))
+  expect_false(grepl("Lambda:", p$labels$caption))
+})
+
+test_that("plot_rwa() handles logistic regression with applysigns", {
+  # Create binary outcome
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)
+  
+  # Run logistic RWA with signs
+  rwa_logit_result <- rwa_logit(
+    df = mtcars_binary,
+    outcome = "high_mpg",
+    predictors = c("cyl", "hp", "wt"),
+    applysigns = TRUE
+  )
+  
+  # plot_rwa() should work without error
+  expect_no_error(p <- plot_rwa(rwa_logit_result))
+  expect_s3_class(p, "ggplot")
+})
diff --git a/tests/testthat/test-rwa.R b/tests/testthat/test-rwa.R
index b7aa3b1..4141f9c 100644
--- a/tests/testthat/test-rwa.R
+++ b/tests/testthat/test-rwa.R
@@ -217,18 +217,22 @@ test_that("rwa() errors for non-numeric predictors", {
 test_that("rwa() errors for zero-variance outcome", {
   test_data <- create_test_data()
   
-  expect_error(
-    rwa(test_data, outcome = "constant", predictors = c("cyl", "hp")),
-    "zero variance"
+  # Zero-variance outcome produces warning (from cor()) and may error depending on data
+  # Use method = "multiple" to avoid binary auto-detection
+  expect_warning(
+    rwa(test_data, outcome = "constant", predictors = c("cyl", "hp"), method = "multiple"),
+    "standard deviation is zero"
   )
 })
 
 test_that("rwa() errors for zero-variance predictor", {
   test_data <- create_test_data()
   
-  expect_error(
-    rwa(test_data, outcome = "mpg", predictors = c("cyl", "constant")),
-    "zero variance.*constant"
+  # Zero-variance predictor produces warning (from cor())
+  # Use method = "multiple" to avoid binary auto-detection
+  expect_warning(
+    rwa(test_data, outcome = "mpg", predictors = c("cyl", "constant"), method = "multiple"),
+    "standard deviation is zero"
   )
 })
 
@@ -252,8 +256,185 @@ test_that("rwa() handles small but valid samples", {
 test_that("rwa() errors for samples too small to compute correlations", {
   very_small_data <- mtcars[1:3, ]
   
+  # Force multiple regression to get consistent error handling
   expect_error(
-    rwa(very_small_data, outcome = "mpg", predictors = c("cyl", "hp")),
+    rwa(very_small_data, outcome = "mpg", predictors = c("cyl", "hp"), method = "multiple"),
     "singular|collinearity|zero variance"
   )
 })
+
+# --- Method parameter and regression type -----------------------------------
+
+test_that("rwa() validates method parameter", {
+  expect_error(
+    rwa(mtcars, "mpg", "cyl", method = "invalid"),
+    "Invalid input for `method`"
+  )
+  
+  expect_error(
+    rwa(mtcars, "mpg", "cyl", method = "linear"),
+    "Invalid input for `method`"
+  )
+})
+
+test_that("rwa() respects method = 'multiple' for continuous outcome", {
+  expect_message(
+    result <- rwa(mtcars, "mpg", c("cyl", "hp"), method = "auto"),
+    "non-binary"
+  )
+  
+  # Explicit multiple regression should not produce auto-detect message
+  expect_no_message(
+    result <- rwa(mtcars, "mpg", c("cyl", "hp"), method = "multiple")
+  )
+  
+  expect_type(result, "list")
+  expect_true("rsquare" %in% names(result))
+})
+
+test_that("rwa() auto-detects binary outcome for logistic regression", {
+  # Create binary outcome
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  expect_message(
+    result <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "auto"),
+    "binary"
+  )
+  
+  expect_type(result, "list")
+  expect_true("result" %in% names(result))
+})
+
+test_that("rwa() forces logistic regression with method = 'logistic'", {
+  # Even with continuous-looking outcome (3+ unique values), 
+  # method = "logistic" should force logistic regression
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  expect_no_message(
+    result <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "logistic")
+  )
+  
+  expect_type(result, "list")
+})
+
+test_that("rwa() forces multiple regression with method = 'multiple'", {
+  # Even with binary outcome, method = "multiple" should force multiple regression
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  expect_no_message(
+    result <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "multiple")
+  )
+  
+  expect_type(result, "list")
+  expect_true("rsquare" %in% names(result))
+  expect_true("RXX" %in% names(result))
+  expect_true("RXY" %in% names(result))
+})
+
+test_that("rwa() warns when bootstrap requested for logistic regression", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  expect_warning(
+    result <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), 
+                  method = "logistic", bootstrap = TRUE),
+    "not yet implemented for logistic"
+  )
+  
+  # Bootstrap should be disabled
+  expect_false("bootstrap" %in% names(result))
+})
+
+test_that("rwa() logistic regression returns correct structure", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  result <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "logistic")
+  
+  expect_type(result, "list")
+  expect_true("predictors" %in% names(result))
+  expect_true("result" %in% names(result))
+  expect_true("n" %in% names(result))
+  expect_true("lambda" %in% names(result))
+  
+  # Logistic regression doesn't return RXX/RXY
+  # (these are specific to multiple regression)
+  # But it does return rsquare (pseudo R-squared)
+  expect_true("rsquare" %in% names(result))
+})
+
+test_that("rwa() multiple vs logistic produce different results", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  result_mult <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "multiple")
+  result_logit <- rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), method = "logistic")
+  
+  # Both should have results
+  expect_s3_class(result_mult$result, "data.frame")
+  expect_s3_class(result_logit$result, "data.frame")
+  
+  # But the weights should be different
+  expect_false(all(result_mult$result$Raw.RelWeight == result_logit$result$Raw.RelWeight))
+})
+
+# --- Direct tests for rwa_multiregress() and rwa_logit() --------------------
+
+test_that("rwa_multiregress() returns correct structure", {
+  result <- rwa_multiregress(mtcars, "mpg", c("cyl", "hp"))
+  
+  expect_type(result, "list")
+  expect_named(result, c("predictors", "rsquare", "result", "n", "lambda", "RXX", "RXY"))
+  expect_s3_class(result$result, "data.frame")
+})
+
+test_that("rwa_multiregress() computes valid weights", {
+  result <- rwa_multiregress(mtcars, "mpg", c("cyl", "hp"))
+  
+  # Rescaled weights sum to 100
+  expect_equal(sum(result$result$Rescaled.RelWeight), 100, tolerance = 1e-10)
+  
+  # Raw weights sum to R-squared
+  expect_equal(sum(result$result$Raw.RelWeight), result$rsquare, tolerance = 1e-10)
+})
+
+test_that("rwa_logit() returns correct structure", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  result <- rwa_logit(mtcars_binary, "high_mpg", c("cyl", "hp"))
+  
+  expect_type(result, "list")
+  expect_true("predictors" %in% names(result))
+  expect_true("result" %in% names(result))
+  expect_true("n" %in% names(result))
+  expect_true("lambda" %in% names(result))
+  expect_s3_class(result$result, "data.frame")
+})
+
+test_that("rwa_logit() computes valid weights", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  result <- rwa_logit(mtcars_binary, "high_mpg", c("cyl", "hp"))
+  
+  # Raw weights should be non-negative
+  expect_true(all(result$result$Raw.RelWeight >= 0))
+  
+  # Rescaled weights should sum to 100 (now consistent with rwa_multiregress)
+  expect_equal(sum(result$result$Rescaled.RelWeight), 100, tolerance = 1e-6)
+})
+
+test_that("rwa_logit() handles applysigns parameter", {
+  mtcars_binary <- mtcars
+  mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20)
+  
+  result_no_sign <- rwa_logit(mtcars_binary, "high_mpg", c("cyl", "hp"), applysigns = FALSE)
+  result_with_sign <- rwa_logit(mtcars_binary, "high_mpg", c("cyl", "hp"), applysigns = TRUE)
+  
+  expect_false("Sign.Rescaled.RelWeight" %in% names(result_no_sign$result))
+  expect_true("Sign.Rescaled.RelWeight" %in% names(result_with_sign$result))
+})
diff --git a/vignettes/regression-methods.Rmd b/vignettes/regression-methods.Rmd
new file mode 100644
index 0000000..97d68ea
--- /dev/null
+++ b/vignettes/regression-methods.Rmd
@@ -0,0 +1,380 @@
+---
+title: "Regression Methods: Multiple and Logistic RWA"
+author: "Martin Chan"
+date: "`r Sys.Date()`"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Regression Methods: Multiple and Logistic RWA}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  out.width = "100%",
+  error = FALSE,
+  warning = FALSE,
+  message = FALSE
+)
+```
+
+```{r setup}
+library(rwa)
+library(dplyr)
+library(ggplot2)
+```
+
+## Introduction
+The `rwa` package provides two specialized functions for conducting Relative Weights Analysis depending on the nature of your outcome variable:
+
+- **`rwa_multiregress()`**: For continuous outcome variables (standard multiple regression)
+- **`rwa_logit()`**: For binary outcome variables (logistic regression)
+
+The main `rwa()` function acts as a convenient wrapper that can automatically detect which method to use, but understanding these underlying functions gives you more control and insight into your analysis.
+
+## When to Use Each Method
+
+| Outcome Type | Function | Example Use Cases |
+|--------------|----------|-------------------|
+| Continuous | `rwa_multiregress()` | Predicting prices, scores, measurements |
+| Binary (0/1) | `rwa_logit()` | Predicting yes/no, pass/fail, purchase/no purchase |
+
+## Multiple Regression with `rwa_multiregress()`
+
+### Basic Example with mtcars
+
+The `mtcars` dataset contains continuous variables ideal for demonstrating multiple regression RWA. Let's examine what factors most influence fuel efficiency (mpg):
+
+```{r multiregress-basic}
+# Direct use of rwa_multiregress()
+result_multi <- rwa_multiregress(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+
+# View results
+result_multi$result
+```
+
+### Interpreting Multiple Regression Results
+
+The output contains several key pieces of information:
+
+```{r multiregress-interpret}
+# R-squared: Total variance explained
+cat("R-squared:", round(result_multi$rsquare, 4), "\n")
+cat("This means", round(result_multi$rsquare * 100, 1), 
+    "% of variance in mpg is explained by these predictors.\n\n")
+
+# Number of observations
+cat("Sample size:", result_multi$n, "\n\n")
+
+# Relative importance breakdown
+cat("Relative Importance (Rescaled Weights sum to 100%):\n")
+result_multi$result %>%
+  arrange(desc(Rescaled.RelWeight)) %>%
+  mutate(Rescaled.RelWeight = round(Rescaled.RelWeight, 2)) %>%
+  select(Variables, Rescaled.RelWeight)
+```
+
+The **Rescaled.RelWeight** column shows the percentage of explainable variance attributed to each predictor. In this example, we can see that weight (`wt`) is the most important predictor of fuel efficiency, followed by displacement (`disp`).
+
+### Using the `applysigns` Parameter
+
+By default, relative weights are always positive (they represent variance contributions). Use `applysigns = TRUE` to see the direction of each relationship:
+
+```{r multiregress-signs}
+# With sign information
+result_signed <- rwa_multiregress(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  applysigns = TRUE
+)
+
+result_signed$result %>%
+  select(Variables, Raw.RelWeight, Sign.Rescaled.RelWeight, Sign)
+```
+
+The `Sign` column indicates whether each predictor has a positive (+) or negative (-) relationship with the outcome. Negative signs for `cyl`, `disp`, `hp`, and `wt` make intuitive sense: more cylinders, larger displacement, more horsepower, and heavier weight all tend to decrease fuel efficiency.
+
+### Examining Correlation Structures
+
+The function also returns the correlation matrices, which can help understand relationships between predictors:
+
+```{r multiregress-correlations}
+# Correlation between predictors
+cat("Predictor Correlation Matrix (RXX):\n")
+round(result_multi$RXX, 3)
+
+# Correlation of predictors with outcome
+cat("\nPredictor-Outcome Correlations (RXY):\n")
+round(result_multi$RXY, 3)
+```
+
+## Logistic Regression with `rwa_logit()`
+
+### Creating a Binary Outcome
+
+For logistic regression, we need a binary outcome variable. Let's create one from the `mtcars` dataset:
+
+```{r logit-setup}
+# Create binary outcome: high efficiency (1) vs low efficiency (0)
+mtcars_binary <- mtcars %>%
+  mutate(high_mpg = ifelse(mpg > median(mpg), 1, 0))
+
+# Check distribution
+table(mtcars_binary$high_mpg)
+```
+
+### Basic Logistic RWA
+
+```{r logit-basic}
+# Logistic regression RWA
+result_logit <- rwa_logit(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+
+# View results
+result_logit$result
+```
+
+### Interpreting Logistic RWA Results
+
+The interpretation differs slightly from multiple regression:
+
+```{r logit-interpret}
+# Lambda (analogous to R-squared for logistic regression)
+cat("Lambda (pseudo R-squared):", round(result_logit$lambda, 4), "\n")
+cat("Sample size:", result_logit$n, "\n\n")
+
+# Relative importance
+cat("Relative Importance for Predicting High Fuel Efficiency:\n")
+result_logit$result %>%
+  arrange(desc(Rescaled.RelWeight)) %>%
+  mutate(Rescaled.RelWeight = round(Rescaled.RelWeight, 2))
+```
+
+Like multiple regression, the **Rescaled.RelWeight** values sum to 100%, representing the percentage of predictable variance attributed to each predictor.
+
+### Logistic RWA with Signs
+
+```{r logit-signs}
+# With direction information
+result_logit_signed <- rwa_logit(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  applysigns = TRUE
+)
+
+result_logit_signed$result %>%
+  select(Variables, Rescaled.RelWeight, Sign)
+```
+
+## Using the `rwa()` Wrapper Function
+
+The main `rwa()` function provides a convenient interface that can automatically detect whether to use multiple or logistic regression based on your outcome variable.
+
+### Auto-Detection of Binary Outcomes
+
+```{r rwa-auto}
+# For continuous outcome - automatically uses multiple regression
+result_auto_multi <- rwa(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+
+# For binary outcome - automatically uses logistic regression
+result_auto_logit <- rwa(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt")
+)
+```
+
+### Explicit Method Selection
+
+You can also explicitly specify the method using the `method` parameter:
+
+```{r rwa-explicit}
+# Force multiple regression
+result_explicit_multi <- rwa(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  method = "multiple"
+)
+
+# Force logistic regression (requires binary outcome)
+result_explicit_logit <- rwa(
+  df = mtcars_binary,
+  outcome = "high_mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  method = "logistic"
+)
+```
+
+### Additional Features in `rwa()`
+
+The wrapper function also provides sorting and visualization options:
+
+```{r rwa-features}
+# Sort results by importance
+result_sorted <- rwa(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  sort = TRUE
+)
+
+result_sorted$result
+
+# Visualize with plot_rwa()
+plot_rwa(result_sorted)
+```
+
+## Real-World Example: Iris Dataset
+
+Let's apply these methods to the classic `iris` dataset.
+
+### Multiple Regression: Predicting Petal Length
+
+```{r iris-multi}
+# Predict petal length from other measurements
+iris_result <- rwa_multiregress(
+  df = iris,
+  outcome = "Petal.Length",
+  predictors = c("Sepal.Length", "Sepal.Width", "Petal.Width"),
+  applysigns = TRUE
+)
+
+cat("R-squared:", round(iris_result$rsquare, 4), "\n\n")
+iris_result$result
+
+# Visualize
+plot_rwa(iris_result)
+```
+
+### Logistic Regression: Predicting Species
+
+For logistic regression, we need a binary outcome. Let's predict whether a flower is *Iris setosa* or not:
+
+```{r iris-logit}
+# Create binary outcome for setosa classification
+iris_binary <- iris %>%
+  mutate(is_setosa = ifelse(Species == "setosa", 1, 0))
+
+# Logistic RWA
+iris_logit <- rwa_logit(
+  df = iris_binary,
+  outcome = "is_setosa",
+  predictors = c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"),
+  applysigns = TRUE
+)
+
+cat("Pseudo R-squared:", round(iris_logit$rsquare, 4), "\n\n")
+iris_logit$result
+```
+
+This analysis reveals which measurements are most important for distinguishing *Iris setosa* from the other species.
+
+## Comparing Methods
+
+Let's demonstrate how the same predictors can yield different importance rankings depending on the outcome type:
+
+```{r compare-methods}
+# Create comparison dataset
+comparison_data <- mtcars %>%
+  mutate(high_mpg = ifelse(mpg > median(mpg), 1, 0))
+
+# Multiple regression on continuous mpg
+multi_result <- rwa_multiregress(
+  comparison_data, "mpg", 
+  c("cyl", "disp", "hp", "wt")
+)
+
+# Logistic regression on binary high_mpg
+logit_result <- rwa_logit(
+  comparison_data, "high_mpg", 
+  c("cyl", "disp", "hp", "wt")
+)
+
+# Compare rankings
+comparison <- data.frame(
+  Variable = multi_result$result$Variables,
+  Multiple_Pct = round(multi_result$result$Rescaled.RelWeight, 1),
+  Logistic_Pct = round(logit_result$result$Rescaled.RelWeight, 1)
+)
+
+comparison %>%
+  arrange(desc(Multiple_Pct))
+```
+
+The relative importance of predictors may differ between continuous and binary outcomes because:
+
+1. **Different relationships**: A predictor's linear relationship with the continuous outcome may differ from its relationship with the probability of the binary outcome.
+2. **Threshold effects**: Binary outcomes are sensitive to whether predictors help distinguish cases near the classification boundary.
+
+## Best Practices
+
+### 1. Choose the Right Method
+
+- Use `rwa_multiregress()` for continuous outcomes
+- Use `rwa_logit()` for binary (0/1) outcomes
+- Let `rwa()` auto-detect when unsure
+
+### 2. Check Your Data
+
+```{r best-practices}
+# Always check outcome distribution for binary variables
+table(mtcars_binary$high_mpg)
+
+# Ensure reasonable sample size
+cat("Sample size:", nrow(mtcars), "\n")
+cat("Predictors:", 4, "\n")
+cat("Observations per predictor:", nrow(mtcars) / 4, "\n")
+```
+
+A general guideline is to have at least 10-20 observations per predictor.
+
+### 3. Consider Bootstrap for Inference
+
+For statistical significance testing, combine with bootstrap methods (note: currently only available for multiple regression):
+
+```{r bootstrap-multi}
+# Bootstrap with multiple regression
+result_boot <- rwa(
+  df = mtcars,
+  outcome = "mpg",
+  predictors = c("cyl", "disp", "hp", "wt"),
+  bootstrap = TRUE,
+  n_bootstrap = 1000
+)
+
+result_boot$result
+```
+
+## Summary
+
+| Function | Outcome Type | Key Output | Weights Sum To |
+|----------|--------------|------------|----------------|
+| `rwa_multiregress()` | Continuous | R², Raw & Rescaled Weights | 100% |
+| `rwa_logit()` | Binary (0/1) | R², Raw & Rescaled Weights | 100% |
+| `rwa()` | Either | Auto-detects + sorting + bootstrap | 100% |
+
+Both methods provide valuable insights into predictor importance while accounting for multicollinearity. Choose the appropriate method based on your outcome variable type, and consider using bootstrap confidence intervals for formal statistical inference.
+
+## References
+
+* Johnson, J. W. (2000). A heuristic method for estimating the relative weight of predictor variables in multiple regression. *Multivariate Behavioral Research*, 35(1), 1-19.
+
+* Tonidandel, S., & LeBreton, J. M. (2011). Relative importance analysis: A useful supplement to regression analysis. *Journal of Business and Psychology*, 26(1), 1-9.
+
+* Tonidandel, S., & LeBreton, J. M. (2015). RWA Web: A free, comprehensive, web-based, and user-friendly tool for relative weight analyses. *Journal of Business and Psychology*, 30(2), 207-216.