diff --git a/NAMESPACE b/NAMESPACE index 29c37a1..6be7c1f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,7 +21,9 @@ importFrom(magrittr,"%>%") importFrom(purrr,map_dfr) importFrom(stats,binomial) importFrom(stats,coef) +importFrom(stats,complete.cases) importFrom(stats,cor) +importFrom(stats,cov.wt) importFrom(stats,glm) importFrom(stats,lm) importFrom(stats,predict) diff --git a/NEWS.md b/NEWS.md index fd6a45f..c3208fd 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,19 @@ # 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. +## New Features + +- Added `rwa_logit()` and `rwa_multiregress()` to support logistic regression and multiple regression. +- Added new vignette to cover the new regression methods. +- Added `use` parameter to `rwa()` function to control how missing data is handled when computing correlations. Options include "pairwise.complete.obs" (default, pairwise deletion), "complete.obs" (listwise deletion), and other standard options from `cor()`. (#12) +- Added `weight` parameter to `rwa()` function to perform weighted Relative Weights Analysis (#12). When a weight variable is specified, the function computes a weighted correlation matrix using `cov.wt()`. This feature allows for proper handling of survey weights or importance weights in the analysis. + +## Improvements + +- Updated all bootstrap functions to support the new `use` and `weight` parameters +- Enhanced input validation to check weight variable properties (numeric, positive) +- Improved test coverage and minor bugfixes + +--- # rwa 0.1.1 diff --git a/R/bootstrap_rwa.R b/R/bootstrap_rwa.R index 1310a6b..3ce3935 100644 --- a/R/bootstrap_rwa.R +++ b/R/bootstrap_rwa.R @@ -24,12 +24,65 @@ NULL #' @param outcome Name of outcome variable #' @param predictors Names of predictor variables #' @param return_all If TRUE, returns list with raw weights, rescaled weights, and rsquare +#' @param use Method for handling missing data in correlations (passed to cor()) +#' @param weight Optional name of weight variable for weighted correlations #' #' @return Numeric vector of raw weights, or list if return_all=TRUE #' @noRd -rwa_core_calculation <- function(thedata, outcome, predictors, return_all = FALSE) { - cor_matrix <- cor(thedata, use = "pairwise.complete.obs") %>% - as.data.frame(stringsAsFactors = FALSE, row.names = NULL) %>% +rwa_core_calculation <- function(thedata, outcome, predictors, return_all = FALSE, + use = "pairwise.complete.obs", weight = NULL) { + + # Compute correlation matrix (weighted or unweighted) + if (!is.null(weight)) { + # Check if weight variable exists in thedata + if (!weight %in% names(thedata)) { + stop(sprintf("Weight variable '%s' not found in data.", weight)) + } + + # Extract weights and variables for analysis + weight_values <- thedata[[weight]] + analysis_data <- thedata %>% dplyr::select(dplyr::all_of(c(outcome, predictors))) + + # Handle NA weights consistently with use parameter + if (any(is.na(weight_values))) { + if (use == "all.obs") { + stop("Weight variable contains NA values and use = 'all.obs'. Set use = 'complete.obs' for listwise deletion.") + } else { + # Remove rows with NA weights + non_na_idx <- !is.na(weight_values) + weight_values <- weight_values[non_na_idx] + analysis_data <- analysis_data[non_na_idx, ] + } + } + + # Require complete cases for cov.wt + complete_cases_idx <- stats::complete.cases(analysis_data) + weight_values <- weight_values[complete_cases_idx] + analysis_data <- analysis_data[complete_cases_idx, , drop = FALSE] + + if (nrow(analysis_data) == 0 || sum(weight_values) == 0) { + stop("Insufficient data for weighted correlation computation.") + } + + # Compute weighted covariance matrix using cov.wt + cov_result <- stats::cov.wt( + x = analysis_data, + wt = weight_values, + cor = TRUE, + method = "unbiased" + ) + + cor_matrix <- cov_result$cor %>% + as.data.frame(stringsAsFactors = FALSE, row.names = NULL) + + } else { + # Unweighted correlation + cor_matrix <- cor(thedata %>% dplyr::select(dplyr::all_of(c(outcome, predictors))), + use = use) %>% + as.data.frame(stringsAsFactors = FALSE, row.names = NULL) + } + + cor_matrix <- cor_matrix %>% remove_all_na_cols() %>% tidyr::drop_na() @@ -78,18 +131,26 @@ rwa_core_calculation <- function(thedata, outcome, predictors, return_all = FALS #' @param indices Bootstrap sample indices (provided by boot::boot) #' @param outcome Outcome variable name #' @param predictors Vector of predictor variable names +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return Numeric vector of raw relative weights #' @keywords internal #' @noRd -rwa_boot_statistic <- function(data, indices, outcome, predictors) { +rwa_boot_statistic <- function(data, indices, outcome, predictors, use = "pairwise.complete.obs", weight_var = NULL) { sample_data <- data[indices, ] - thedata <- sample_data %>% - dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% - tidyr::drop_na(dplyr::all_of(outcome)) + if (!is.null(weight_var)) { + thedata <- sample_data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors, weight_var))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } else { + thedata <- sample_data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } - rwa_core_calculation(thedata, outcome, predictors, return_all = FALSE) + rwa_core_calculation(thedata, outcome, predictors, return_all = FALSE, use = use, weight = weight_var) } #' Bootstrap statistic function for rescaled RWA weights @@ -100,18 +161,26 @@ rwa_boot_statistic <- function(data, indices, outcome, predictors) { #' @param indices Bootstrap sample indices (provided by boot::boot) #' @param outcome Outcome variable name #' @param predictors Vector of predictor variable names +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return Numeric vector of rescaled relative weights (summing to 100) #' @keywords internal #' @noRd -rwa_boot_statistic_rescaled <- function(data, indices, outcome, predictors) { +rwa_boot_statistic_rescaled <- function(data, indices, outcome, predictors, use = "pairwise.complete.obs", weight_var = NULL) { sample_data <- data[indices, ] - thedata <- sample_data %>% - dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% - tidyr::drop_na(dplyr::all_of(outcome)) + if (!is.null(weight_var)) { + thedata <- sample_data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors, weight_var))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } else { + thedata <- sample_data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } - result <- rwa_core_calculation(thedata, outcome, predictors, return_all = TRUE) + result <- rwa_core_calculation(thedata, outcome, predictors, return_all = TRUE, use = use, weight = weight_var) result$rescaled_weights } @@ -125,23 +194,25 @@ rwa_boot_statistic_rescaled <- function(data, indices, outcome, predictors) { #' @param outcome Outcome variable name #' @param predictors Vector of predictor variable names #' @param focal Focal variable for comparisons (optional) +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return Numeric vector containing raw weights, random comparison differences, #' and (if focal specified) focal comparison differences #' @keywords internal #' @noRd -rwa_boot_comprehensive <- function(data, indices, outcome, predictors, focal = NULL) { +rwa_boot_comprehensive <- function(data, indices, outcome, predictors, focal = NULL, use = "pairwise.complete.obs", weight_var = NULL) { sample_data <- data[indices, ] # Get raw weights using core calculation - raw_weights <- rwa_boot_statistic(sample_data, seq_len(nrow(sample_data)), outcome, predictors) + raw_weights <- rwa_boot_statistic(sample_data, seq_len(nrow(sample_data)), outcome, predictors, use = use, weight_var = weight_var) # Get random variable comparison (difference from random variable) - rand_diff <- rwa_rand_internal(sample_data, outcome, predictors) + rand_diff <- rwa_rand_internal(sample_data, outcome, predictors, use = use, weight = weight_var) # Get focal variable comparison if focal is specified if (!is.null(focal)) { - focal_diff <- rwa_comp_internal(sample_data, outcome, predictors, focal) + focal_diff <- rwa_comp_internal(sample_data, outcome, predictors, focal, use = use, weight = weight_var) c(raw_weights, rand_diff, focal_diff) } else { c(raw_weights, rand_diff) @@ -157,19 +228,28 @@ rwa_boot_comprehensive <- function(data, indices, outcome, predictors, focal = N #' @param df Data frame #' @param outcome Outcome variable name #' @param predictors Vector of predictor variable names +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return Numeric vector of weight differences (predictor weight - random weight) #' @keywords internal #' @noRd -rwa_rand_internal <- function(df, outcome, predictors) { - thedata <- df %>% - dplyr::select(all_of(c(outcome, predictors))) %>% - tidyr::drop_na(all_of(outcome)) %>% - dplyr::mutate(rand = rnorm(dplyr::n(), 0, 1)) +rwa_rand_internal <- function(df, outcome, predictors, use = "pairwise.complete.obs", weight = NULL) { + if (!is.null(weight)) { + thedata <- df %>% + dplyr::select(all_of(c(outcome, predictors, weight))) %>% + tidyr::drop_na(all_of(outcome)) %>% + dplyr::mutate(rand = rnorm(dplyr::n(), 0, 1)) + } else { + thedata <- df %>% + dplyr::select(all_of(c(outcome, predictors))) %>% + tidyr::drop_na(all_of(outcome)) %>% + dplyr::mutate(rand = rnorm(dplyr::n(), 0, 1)) + } # Use core calculation with random variable added predictors_with_rand <- c(predictors, "rand") - RawWgt <- rwa_core_calculation(thedata, outcome, predictors_with_rand, return_all = FALSE) + RawWgt <- rwa_core_calculation(thedata, outcome, predictors_with_rand, return_all = FALSE, use = use, weight = weight) RawWgt <- RawWgt - tail(RawWgt, n = 1) # subtract random variable weight head(RawWgt, -1) # remove random variable from output @@ -184,19 +264,28 @@ rwa_rand_internal <- function(df, outcome, predictors) { #' @param outcome Outcome variable name #' @param predictors Vector of predictor variable names #' @param focal Name of the focal variable to compare against +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return Numeric vector of weight differences (predictor weight - focal weight) #' @keywords internal #' @noRd -rwa_comp_internal <- function(df, outcome, predictors, focal) { - thedata <- df %>% - dplyr::select(all_of(c(outcome, predictors))) %>% - tidyr::drop_na(all_of(outcome)) %>% - dplyr::relocate(all_of(focal), .after = dplyr::last_col()) +rwa_comp_internal <- function(df, outcome, predictors, focal, use = "pairwise.complete.obs", weight = NULL) { + if (!is.null(weight)) { + thedata <- df %>% + dplyr::select(all_of(c(outcome, predictors, weight))) %>% + tidyr::drop_na(all_of(outcome)) %>% + dplyr::relocate(all_of(focal), .after = dplyr::last_col()) + } else { + thedata <- df %>% + dplyr::select(all_of(c(outcome, predictors))) %>% + tidyr::drop_na(all_of(outcome)) %>% + dplyr::relocate(all_of(focal), .after = dplyr::last_col()) + } # Reorder predictors to match data predictors_reordered <- c(predictors[predictors != focal], focal) - RawWgt <- rwa_core_calculation(thedata, outcome, predictors_reordered, return_all = FALSE) + RawWgt <- rwa_core_calculation(thedata, outcome, predictors_reordered, return_all = FALSE, use = use, weight = weight) RawWgt <- RawWgt - tail(RawWgt, n = 1) # subtract focal variable weight head(RawWgt, -1) # remove focal variable from output @@ -304,6 +393,8 @@ extract_ci <- function(boot_object, conf_level = 0.95, variable_names = NULL, ci #' @param comprehensive Whether to run comprehensive analysis with random #' variable and focal comparisons #' @param include_rescaled Whether to bootstrap rescaled weights +#' @param use Method for handling missing data in correlations +#' @param weight Optional name of weight variable #' #' @return List containing: #' - boot_object: Raw bootstrap object @@ -313,12 +404,18 @@ extract_ci <- function(boot_object, conf_level = 0.95, variable_names = NULL, ci #' @noRd run_rwa_bootstrap <- function(data, outcome, predictors, n_bootstrap = 1000, conf_level = 0.95, focal = NULL, comprehensive = FALSE, - include_rescaled = FALSE) { + include_rescaled = FALSE, use = "pairwise.complete.obs", weight = NULL) { # Prepare data - bootstrap_data <- data %>% - dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% - tidyr::drop_na(dplyr::all_of(outcome)) + if (!is.null(weight)) { + bootstrap_data <- data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors, weight))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } else { + bootstrap_data <- data %>% + dplyr::select(dplyr::all_of(c(outcome, predictors))) %>% + tidyr::drop_na(dplyr::all_of(outcome)) + } # Check sample size if (nrow(bootstrap_data) < 50) { @@ -326,12 +423,15 @@ run_rwa_bootstrap <- function(data, outcome, predictors, n_bootstrap = 1000, } # Always bootstrap raw weights + # Note: using weight_var to avoid partial matching with boot::boot's "weights" parameter boot_result_raw <- boot::boot( data = bootstrap_data, statistic = rwa_boot_statistic, R = n_bootstrap, outcome = outcome, - predictors = predictors + predictors = predictors, + use = use, + weight_var = weight ) # Extract CIs for raw weights @@ -343,12 +443,15 @@ run_rwa_bootstrap <- function(data, outcome, predictors, n_bootstrap = 1000, # Bootstrap rescaled weights if requested if (include_rescaled) { + # Note: using weight_var to avoid partial matching with boot::boot's "weights" parameter boot_result_rescaled <- boot::boot( data = bootstrap_data, statistic = rwa_boot_statistic_rescaled, R = n_bootstrap, outcome = outcome, - predictors = predictors + predictors = predictors, + use = use, + weight_var = weight ) rescaled_ci <- extract_ci(boot_result_rescaled, conf_level, predictors, "rescaled") @@ -358,13 +461,16 @@ run_rwa_bootstrap <- function(data, outcome, predictors, n_bootstrap = 1000, # Handle comprehensive analysis if requested if (comprehensive && !is.null(focal)) { + # Note: using weight_var to avoid partial matching with boot::boot's "weights" parameter boot_result_comp <- boot::boot( data = bootstrap_data, statistic = rwa_boot_comprehensive, R = n_bootstrap, outcome = outcome, predictors = predictors, - focal = focal + focal = focal, + use = use, + weight_var = weight ) n_vars <- length(predictors) diff --git a/R/rwa.R b/R/rwa.R index 05c0f7a..25c37b7 100644 --- a/R/rwa.R +++ b/R/rwa.R @@ -46,6 +46,19 @@ #' @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 use Method for handling missing data when computing correlations. Options are: +#' "everything" (missing values in correlations propagate), +#' "all.obs" (error if missing values present), +#' "complete.obs" (listwise deletion), +#' "na.or.complete" (error if some but not all missing), +#' "pairwise.complete.obs" (pairwise deletion, default). +#' See \code{\link[stats]{cor}} for more details. Only applicable for multiple regression. +#' Note: When \code{weight} is specified, complete cases (listwise deletion) is +#' always used for weighted correlation computation regardless of \code{use}. +#' @param weight Optional name of a weight variable in the data frame. If provided, +#' a weighted correlation matrix will be computed using the specified weights. +#' The weight variable must be numeric and positive. Defaults to \code{NULL} +#' (unweighted analysis). Only applicable for multiple regression. #' #' @return `rwa()` returns a list of outputs, as follows: #' - `predictors`: character vector of names of the predictor variables used. @@ -93,6 +106,14 @@ #' # For faster examples, use a subset of data for bootstrap #' diamonds_small <- diamonds[sample(nrow(diamonds), 1000), ] #' +#' # RWA with different missing data handling +#' # Use complete.obs for listwise deletion +#' rwa(diamonds_small, "price", c("depth", "carat"), use = "complete.obs") +#' +#' # RWA with weights +#' diamonds_small$sample_weight <- runif(nrow(diamonds_small), 0.5, 2) +#' rwa(diamonds_small, "price", c("depth", "carat"), weight = "sample_weight") +#' #' # RWA with bootstrap confidence intervals (raw weights only) #' rwa(diamonds_small, "price", c("depth", "carat"), #' bootstrap = TRUE, n_bootstrap = 100) @@ -125,8 +146,9 @@ rwa <- function(df, conf_level = 0.95, focal = NULL, comprehensive = FALSE, - include_rescaled_ci = FALSE) { - + include_rescaled_ci = FALSE, + use = "pairwise.complete.obs", + weight = NULL) { # ---- Input validation ---- @@ -148,6 +170,30 @@ rwa <- function(df, stop("`n_bootstrap` must be a positive integer.") } + # Validate use parameter + valid_use_options <- c("everything", "all.obs", "complete.obs", + "na.or.complete", "pairwise.complete.obs") + if (!use %in% valid_use_options) { + stop(sprintf("`use` must be one of: %s", + paste(valid_use_options, collapse = ", "))) + } + + # Validate weight parameter if provided + if (!is.null(weight)) { + if (!is.character(weight) || length(weight) != 1) { + stop("`weight` must be a single character string specifying the weight variable name.") + } + if (!weight %in% names(df)) { + stop(sprintf("Weight variable '%s' not found in data.", weight)) + } + if (!is.numeric(df[[weight]])) { + stop(sprintf("Weight variable '%s' must be numeric.", weight)) + } + if (any(df[[weight]] <= 0, na.rm = TRUE)) { + stop(sprintf("Weight variable '%s' must have positive values.", weight)) + } + } + # Check that outcome and predictors exist in data if (!outcome %in% names(df)) { stop(sprintf("Outcome variable '%s' not found in data.", outcome)) @@ -201,6 +247,13 @@ rwa <- function(df, bootstrap <- FALSE } + # ---- Handle weight/use parameters for logistic regression ---- + + if (use_logistic && (!is.null(weight) || use != "pairwise.complete.obs")) { + warning("Weight and use parameters are only applicable for multiple regression. ", + "They will be ignored for logistic regression.") + } + # ---- Call appropriate sub-function ---- if (use_logistic) { @@ -215,7 +268,9 @@ rwa <- function(df, df = df, outcome = outcome, predictors = predictors, - applysigns = applysigns + applysigns = applysigns, + use = use, + weight = weight ) } @@ -239,7 +294,9 @@ rwa <- function(df, conf_level = conf_level, focal = focal, comprehensive = comprehensive, - include_rescaled = include_rescaled_ci + include_rescaled = include_rescaled_ci, + use = use, + weight = weight ) # Add confidence intervals to result dataframe diff --git a/R/rwa_multiregress.R b/R/rwa_multiregress.R index 4f2f108..c9e144d 100644 --- a/R/rwa_multiregress.R +++ b/R/rwa_multiregress.R @@ -10,12 +10,25 @@ #' `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. +#' See 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`. +#' @param use Method for handling missing data when computing correlations. Options are: +#' "everything" (missing values in correlations propagate), +#' "all.obs" (error if missing values present), +#' "complete.obs" (listwise deletion), +#' "na.or.complete" (error if some but not all missing), +#' "pairwise.complete.obs" (pairwise deletion, default). +#' See \code{\link[stats]{cor}} for more details. +#' Note: When \code{weight} is specified, complete cases (listwise deletion) is +#' always used for weighted correlation computation regardless of \code{use}. +#' @param weight Optional name of a weight variable in the data frame. If provided, +#' a weighted correlation matrix will be computed using the specified weights. +#' The weight variable must be numeric and positive. Defaults to \code{NULL} +#' (unweighted analysis). #' #' @return `rwa_multiregress()` returns a list of outputs, as follows: #' - `predictors`: character vector of names of the predictor variables used. @@ -30,7 +43,7 @@ #' #' @importFrom magrittr %>% #' @importFrom tidyr drop_na -#' @importFrom stats cor +#' @importFrom stats cor cov.wt complete.cases #' @import dplyr #' @examples #' # Basic multiple regression RWA @@ -52,21 +65,101 @@ #' ) #' result_signed$result #' +#' # Using listwise deletion for missing data +#' rwa_multiregress( +#' df = mtcars, +#' outcome = "mpg", +#' predictors = c("cyl", "disp"), +#' use = "complete.obs" +#' ) +#' +#' # With observation weights +#' mtcars_weighted <- mtcars +#' mtcars_weighted$w <- runif(nrow(mtcars), 0.5, 2) +#' rwa_multiregress( +#' df = mtcars_weighted, +#' outcome = "mpg", +#' predictors = c("cyl", "disp"), +#' weight = "w" +#' ) +#' #' @export rwa_multiregress <- function(df, outcome, predictors, - applysigns = FALSE){ + applysigns = FALSE, + use = "pairwise.complete.obs", + weight = NULL){ # Gets data frame in right order and form - thedata <- - df %>% - dplyr::select(all_of(c(outcome, predictors))) %>% - tidyr::drop_na(all_of(outcome)) + if (!is.null(weight)) { + thedata <- + df %>% + dplyr::select(all_of(c(outcome, predictors, weight))) %>% + tidyr::drop_na(all_of(outcome)) + } else { + 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) %>% + # Compute correlation matrix (weighted or unweighted) + if (!is.null(weight)) { + # Extract weights and variables for analysis + weight_values <- thedata[[weight]] + analysis_data <- thedata %>% dplyr::select(all_of(c(outcome, predictors))) + + # Handle NA weights consistently with use parameter + if (any(is.na(weight_values))) { + if (use == "all.obs") { + stop("Weight variable contains NA values and use = 'all.obs'. Set use = 'complete.obs' for listwise deletion.") + } else { + # Remove rows with NA weights + non_na_idx <- !is.na(weight_values) + weight_values <- weight_values[non_na_idx] + analysis_data <- analysis_data[non_na_idx, ] + } + } + + if (sum(weight_values) == 0) { + stop("Sum of weights is zero. Cannot compute weighted correlation.") + } + + # Compute weighted covariance matrix using cov.wt + # Note: cov.wt requires complete cases + complete_cases_idx <- stats::complete.cases(analysis_data) + if (sum(complete_cases_idx) == 0) { + stop("No complete cases available for weighted correlation computation.") + } + + # Track actual n used for weighted analysis + n_used <- sum(complete_cases_idx) + + cov_result <- stats::cov.wt( + x = analysis_data[complete_cases_idx, , drop = FALSE], + wt = weight_values[complete_cases_idx], + cor = TRUE, + method = "unbiased" + ) + + cor_matrix <- cov_result$cor %>% + as.data.frame(stringsAsFactors = FALSE, row.names = NULL) + + } else { + # Unweighted correlation + cor_matrix <- + stats::cor(thedata[, c(outcome, predictors)], use = use) %>% + as.data.frame(stringsAsFactors = FALSE, row.names = NULL) + + # Track n for unweighted analysis (complete cases on all variables). + # Note: When use = "pairwise.complete.obs", individual correlations may + # use more observations than reported here. n reflects the most + # conservative count (complete cases across all variables). + n_used <- nrow(tidyr::drop_na(thedata)) + } + + cor_matrix <- cor_matrix %>% remove_all_na_cols() %>% tidyr::drop_na() @@ -111,8 +204,6 @@ rwa_multiregress <- function(df, Rescaled.RelWeight = import, Sign = sign) # Output - results - complete_cases <- nrow(tidyr::drop_na(thedata)) - if(applysigns == TRUE){ result <- result %>% @@ -124,7 +215,7 @@ rwa_multiregress <- function(df, list("predictors" = Variables, "rsquare" = rsquare, "result" = result, - "n" = complete_cases, + "n" = n_used, "lambda" = lambda, "RXX" = RXX, "RXY" = RXY) diff --git a/man/rwa.Rd b/man/rwa.Rd index a01aa8c..eddb635 100644 --- a/man/rwa.Rd +++ b/man/rwa.Rd @@ -16,7 +16,9 @@ rwa( conf_level = 0.95, focal = NULL, comprehensive = FALSE, - include_rescaled_ci = FALSE + include_rescaled_ci = FALSE, + use = "pairwise.complete.obs", + weight = NULL ) } \arguments{ @@ -60,6 +62,21 @@ 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{use}{Method for handling missing data when computing correlations. Options are: +"everything" (missing values in correlations propagate), +"all.obs" (error if missing values present), +"complete.obs" (listwise deletion), +"na.or.complete" (error if some but not all missing), +"pairwise.complete.obs" (pairwise deletion, default). +See \code{\link[stats]{cor}} for more details. Only applicable for multiple regression. +Note: When \code{weight} is specified, complete cases (listwise deletion) is +always used for weighted correlation computation regardless of \code{use}.} + +\item{weight}{Optional name of a weight variable in the data frame. If provided, +a weighted correlation matrix will be computed using the specified weights. +The weight variable must be numeric and positive. Defaults to \code{NULL} +(unweighted analysis). Only applicable for multiple regression.} } \value{ \code{rwa()} returns a list of outputs, as follows: @@ -127,6 +144,14 @@ diamonds |> # For faster examples, use a subset of data for bootstrap diamonds_small <- diamonds[sample(nrow(diamonds), 1000), ] +# RWA with different missing data handling +# Use complete.obs for listwise deletion +rwa(diamonds_small, "price", c("depth", "carat"), use = "complete.obs") + +# RWA with weights +diamonds_small$sample_weight <- runif(nrow(diamonds_small), 0.5, 2) +rwa(diamonds_small, "price", c("depth", "carat"), weight = "sample_weight") + # RWA with bootstrap confidence intervals (raw weights only) rwa(diamonds_small, "price", c("depth", "carat"), bootstrap = TRUE, n_bootstrap = 100) diff --git a/man/rwa_multiregress.Rd b/man/rwa_multiregress.Rd index da3353a..b547d20 100644 --- a/man/rwa_multiregress.Rd +++ b/man/rwa_multiregress.Rd @@ -4,7 +4,14 @@ \alias{rwa_multiregress} \title{Create a Relative Weights Analysis (RWA)} \usage{ -rwa_multiregress(df, outcome, predictors, applysigns = FALSE) +rwa_multiregress( + df, + outcome, + predictors, + applysigns = FALSE, + use = "pairwise.complete.obs", + weight = NULL +) } \arguments{ \item{df}{Data frame or tibble to be passed through.} @@ -14,6 +21,21 @@ rwa_multiregress(df, outcome, predictors, applysigns = FALSE) \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{use}{Method for handling missing data when computing correlations. Options are: +"everything" (missing values in correlations propagate), +"all.obs" (error if missing values present), +"complete.obs" (listwise deletion), +"na.or.complete" (error if some but not all missing), +"pairwise.complete.obs" (pairwise deletion, default). +See \code{\link[stats]{cor}} for more details. +Note: When \code{weight} is specified, complete cases (listwise deletion) is +always used for weighted correlation computation regardless of \code{use}.} + +\item{weight}{Optional name of a weight variable in the data frame. If provided, +a weighted correlation matrix will be computed using the specified weights. +The weight variable must be numeric and positive. Defaults to \code{NULL} +(unweighted analysis).} } \value{ \code{rwa_multiregress()} returns a list of outputs, as follows: @@ -42,7 +64,7 @@ maximally related to the original set of predictors. \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. +See \url{https://www.scotttonidandel.com/rwa-web} for the original implementation that inspired this package. } \examples{ # Basic multiple regression RWA @@ -64,4 +86,22 @@ result_signed <- rwa_multiregress( ) result_signed$result +# Using listwise deletion for missing data +rwa_multiregress( + df = mtcars, + outcome = "mpg", + predictors = c("cyl", "disp"), + use = "complete.obs" +) + +# With observation weights +mtcars_weighted <- mtcars +mtcars_weighted$w <- runif(nrow(mtcars), 0.5, 2) +rwa_multiregress( + df = mtcars_weighted, + outcome = "mpg", + predictors = c("cyl", "disp"), + weight = "w" +) + } diff --git a/tests/testthat/test-rwa.R b/tests/testthat/test-rwa.R index 4141f9c..0e658bd 100644 --- a/tests/testthat/test-rwa.R +++ b/tests/testthat/test-rwa.R @@ -438,3 +438,218 @@ test_that("rwa_logit() handles applysigns parameter", { expect_false("Sign.Rescaled.RelWeight" %in% names(result_no_sign$result)) expect_true("Sign.Rescaled.RelWeight" %in% names(result_with_sign$result)) }) + +# --- Missing data handling options ------------------------------------------ + +test_that("rwa() accepts different use parameter values", { + # Test with complete.obs (listwise deletion) + result_complete <- rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), + use = "complete.obs", method = "multiple") + expect_type(result_complete, "list") + expect_equal(result_complete$n, nrow(mtcars)) + + # Test with pairwise.complete.obs (default) + result_pairwise <- rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), + use = "pairwise.complete.obs", method = "multiple") + expect_type(result_pairwise, "list") + expect_equal(result_pairwise$n, nrow(mtcars)) +}) + +test_that("rwa() validates use parameter", { + expect_error( + rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), use = "invalid"), + "use.*must be one of" + ) +}) + +test_that("rwa() handles pairwise vs complete deletion differently with missing data", { + # Create data with missing values in predictors + mtcars_na <- mtcars + mtcars_na$cyl[1:2] <- NA + mtcars_na$hp[3:4] <- NA + + # With pairwise deletion, it should use all available pairwise correlations + result_pairwise <- rwa(mtcars_na, outcome = "mpg", predictors = c("cyl", "hp"), + use = "pairwise.complete.obs", method = "multiple") + + # With complete.obs, it should only use rows with no missing values + result_complete <- rwa(mtcars_na, outcome = "mpg", predictors = c("cyl", "hp"), + use = "complete.obs", method = "multiple") + + # Both should return valid results + expect_type(result_pairwise, "list") + expect_type(result_complete, "list") + + # Both methods use listwise deletion on the outcome variable. + # They differ in how predictor missingness is handled during correlation computation. + # n reports complete cases across all variables (conservative estimate for pairwise). + expect_true(result_pairwise$n <= nrow(mtcars_na)) + expect_true(result_complete$n <= result_pairwise$n) +}) + +# --- Weight variable support ------------------------------------------------ + +test_that("rwa() accepts weight parameter", { + # Add a weight variable + mtcars_weighted <- mtcars + mtcars_weighted$weights <- runif(nrow(mtcars), 0.5, 2) + + result <- rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), + weight = "weights", method = "multiple") + + expect_type(result, "list") + expect_named(result, c("predictors", "rsquare", "result", "n", "lambda", "RXX", "RXY")) + + # Rescaled weights should still sum to 100 + expect_equal(sum(result$result$Rescaled.RelWeight), 100, tolerance = 1e-10) +}) + +test_that("rwa() validates weight parameter", { + mtcars_weighted <- mtcars + mtcars_weighted$weights <- runif(nrow(mtcars), 0.5, 2) + mtcars_weighted$char_weight <- as.character(mtcars_weighted$weights) + mtcars_weighted$neg_weight <- -1 * mtcars_weighted$weights + mtcars_weighted$zero_weight <- rep(0, nrow(mtcars)) + + # Non-existent weight variable + expect_error( + rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), weight = "nonexistent"), + "Weight variable.*not found" + ) + + # Non-numeric weight variable + expect_error( + rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), weight = "char_weight"), + "Weight variable.*must be numeric" + ) + + # Negative weight values + expect_error( + rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), weight = "neg_weight"), + "Weight variable.*must have positive values" + ) + + # Zero weight values (should also be rejected as they're not positive) + expect_error( + rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), weight = "zero_weight"), + "Weight variable.*must have positive values" + ) +}) + +test_that("rwa() produces different results with and without weights", { + # Create data with weights that emphasize certain observations + mtcars_weighted <- mtcars + mtcars_weighted$weights <- rep(1, nrow(mtcars)) + mtcars_weighted$weights[1:10] <- 5 # Give more weight to first 10 observations + + result_unweighted <- rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), method = "multiple") + result_weighted <- rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), + weight = "weights", method = "multiple") + + # Both should return valid results + expect_type(result_unweighted, "list") + expect_type(result_weighted, "list") + + # Results should be different (weights should affect the analysis) + expect_false(isTRUE(all.equal(result_unweighted$result$Raw.RelWeight, + result_weighted$result$Raw.RelWeight))) +}) + +test_that("rwa() handles weights with equal values (equivalent to unweighted)", { + # All equal weights should give same result as unweighted + mtcars_weighted <- mtcars + mtcars_weighted$weights <- rep(1, nrow(mtcars)) + + result_unweighted <- rwa(mtcars, outcome = "mpg", predictors = c("cyl", "hp"), method = "multiple") + result_weighted <- rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), + weight = "weights", method = "multiple") + + # Results should be very similar (allowing for numerical precision) + expect_equal(result_unweighted$result$Raw.RelWeight, + result_weighted$result$Raw.RelWeight, + tolerance = 1e-6) +}) + +test_that("rwa() warns when weight/use used with logistic regression", { + mtcars_binary <- mtcars + mtcars_binary$high_mpg <- as.numeric(mtcars_binary$mpg > 20) + mtcars_binary$weights <- runif(nrow(mtcars), 0.5, 2) + + # Weight parameter with logistic regression should warn + expect_warning( + rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), weight = "weights", method = "logistic"), + "Weight and use parameters are only applicable for multiple regression" + ) + + # use parameter with logistic regression should warn + expect_warning( + rwa(mtcars_binary, "high_mpg", c("cyl", "hp"), use = "complete.obs", method = "logistic"), + "Weight and use parameters are only applicable for multiple regression" + ) +}) + +# --- Parameter combination tests -------------------------------------------- + +test_that("rwa() works with bootstrap and weight combined", { + skip_on_cran() + + mtcars_weighted <- mtcars + mtcars_weighted$weights <- runif(nrow(mtcars), 0.5, 2) + + result <- rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp"), + weight = "weights", bootstrap = TRUE, n_bootstrap = 50, + method = "multiple") + + # Should have bootstrap results + expect_true("bootstrap" %in% names(result)) + expect_true("Raw.RelWeight.CI.Lower" %in% names(result$result)) + expect_true("Raw.RelWeight.CI.Upper" %in% names(result$result)) + + # Results should be valid + expect_equal(sum(result$result$Rescaled.RelWeight), 100, tolerance = 1e-10) +}) + +test_that("rwa() works with bootstrap and use combined", { + skip_on_cran() + + # Create data with some NA values + mtcars_na <- mtcars + mtcars_na$cyl[1:2] <- NA + + result <- rwa(mtcars_na, outcome = "mpg", predictors = c("cyl", "hp"), + use = "complete.obs", bootstrap = TRUE, n_bootstrap = 50, + method = "multiple") + + # Should have bootstrap results + expect_true("bootstrap" %in% names(result)) + expect_true("Raw.RelWeight.CI.Lower" %in% names(result$result)) + + # Sample size should reflect listwise deletion + expect_equal(result$n, nrow(mtcars) - 2) +}) + +test_that("rwa() works with all multiple regression parameters combined", { + skip_on_cran() + + mtcars_weighted <- mtcars + mtcars_weighted$weights <- runif(nrow(mtcars), 0.5, 2) + + # Test with weight + use + bootstrap + sort + applysigns + result <- rwa(mtcars_weighted, outcome = "mpg", predictors = c("cyl", "hp", "wt"), + weight = "weights", + use = "complete.obs", + bootstrap = TRUE, + n_bootstrap = 50, + sort = TRUE, + applysigns = TRUE, + method = "multiple") + + # Should have all expected components + expect_true("bootstrap" %in% names(result)) + expect_true("Sign.Rescaled.RelWeight" %in% names(result$result)) + expect_true("Raw.RelWeight.CI.Lower" %in% names(result$result)) + + # Results should be sorted by Rescaled.RelWeight in descending order + weights <- result$result$Rescaled.RelWeight + expect_true(all(diff(weights) <= 0)) +})