diff --git a/.Rbuildignore b/.Rbuildignore index d013b93..bd119d4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,6 @@ ^_pkgdown\.yml$ ^docs$ ^pkgdown$ +^.*\.Rhistory$ +^_scratch$ +^Rplots\.pdf$ diff --git a/.gitignore b/.gitignore index 221acc4..39b4def 100644 --- a/.gitignore +++ b/.gitignore @@ -3,7 +3,11 @@ inst/doc/ *.Rcheck/ *.tar.gz *.md +!README.md +!NEWS.md +!cran-comments.md figure/ +Rplots.pdf # R user files .Rhistory diff --git a/DESCRIPTION b/DESCRIPTION index a5b61eb..1796b0d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: svyTable1 Title: Create Survey-Weighted Descriptive Statistics and Diagnostic Tables -Version: 0.14.1 +Version: 0.14.1.9000 Authors@R: c(person("Ehsan", "Karim", email = "ehsan.karim@gmail.com", role = c("aut", "cre")), @@ -10,20 +10,28 @@ Authors@R: c(person("Ehsan", "Karim", person("Kevin", "Hu", role = "ctb", comment = "Reported a scoping bug in the svycoxph_CE function regarding survey design objects created via pipes.")) -Description: A tool to create publication-ready tables from complex - survey data, including descriptive summaries and multi-panel - interaction reports. It also provides a suite of functions to - evaluate survey-weighted regression models and survey-weighted - survival plots, including coefficient diagnostics, - goodness-of-fit tests, design-correct AUC calculations, - proportional hazards assumption testing, and visualization - of interaction effects. +Description: A teaching toolkit for graduate epidemiology students analyzing + complex survey data such as 'NHANES'. It builds survey-weighted "Table 1" + summaries with reliability suppression, fits and diagnoses survey-weighted + regression and survival models, computes additive interaction measures, and + pools multiply-imputed analyses. Methods are described in Lumley (2004) + for survey-weighted estimation, Archer and + Lemeshow (2006) for survey + goodness-of-fit, Parker et al. (2017) + for the National + Center for Health Statistics data-presentation reliability standards, + VanderWeele and Knol (2014) and Hosmer and + Lemeshow (1992) for additive + interaction, van Buuren and Groothuis-Oudshoorn (2011) + and Rubin (1987, ISBN:9780471087052) for + multiple imputation, and Westreich and Greenland (2013) + for interpreting adjusted estimates. License: MIT + file LICENSE -URL: https://ehsanx.github.io/svyTable1, https://github.com/ehsanx/svyTable1 +URL: https://ehsanx.github.io/svyTable1/, https://github.com/ehsanx/svyTable1 BugReports: https://github.com/ehsanx/svyTable1/issues Encoding: UTF-8 RoxygenNote: 7.3.3 -Imports: +Imports: dplyr, emmeans, ggplot2, @@ -48,12 +56,12 @@ Imports: utils, WeightedROC, broom -Suggests: +Suggests: NHANES, mitools, testthat (>= 3.0.0) VignetteBuilder: knitr Config/testthat/edition: 3 -Depends: +Depends: R (>= 3.5) LazyData: true diff --git a/NAMESPACE b/NAMESPACE index 4780458..a46204f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -88,7 +88,9 @@ importFrom(stats,model.matrix) importFrom(stats,na.omit) importFrom(stats,plogis) importFrom(stats,pnorm) +importFrom(stats,qlogis) importFrom(stats,qnorm) +importFrom(stats,qt) importFrom(stats,quantile) importFrom(stats,residuals) importFrom(stats,terms) diff --git a/R/.Rhistory b/R/.Rhistory deleted file mode 100644 index b1fc37a..0000000 --- a/R/.Rhistory +++ /dev/null @@ -1,512 +0,0 @@ -addintlist <- function(model, -factor1_name, -factor2_name, -measures = "all", -conf.level = 0.95) { -# Check if the core calculation function 'addint' exists -if (!exists("addint") || !is.function(addint)) { -stop("The 'addint' function is required but was not found. Ensure it is loaded.", call. = FALSE) -} -if (!inherits(model, c("svycoxph", "svyglm", "coxph", "glm"))) { -warning("Model object may not be of a supported type.") -# Attempt to proceed anyway -} -beta <- tryCatch(stats::coef(model), error = function(e) NULL) -if (is.null(beta)) { -warning("Could not extract coefficients from the model.") -return(NULL) -} -coef_names_all <- names(beta) -# --- Get Factor Levels (excluding reference) --- -levels1 <- NULL -levels2 <- NULL -# Try accessing xlevels if they exist (common in lm, glm) -if (!is.null(model$xlevels)) { -levels1 <- model$xlevels[[factor1_name]] -levels2 <- model$xlevels[[factor2_name]] -} else if (inherits(model, "survey.design")) { -# For survey design objects passed directly (less common) -if (factor1_name %in% names(model$variables) && is.factor(model$variables[[factor1_name]])) { -levels1 <- levels(model$variables[[factor1_name]]) -} -if (factor2_name %in% names(model$variables) && is.factor(model$variables[[factor2_name]])) { -levels2 <- levels(model$variables[[factor2_name]]) -} -} else if (!is.null(model$survey.design)) { -# For survey model objects (svyglm, svycoxph) -model_data_vars <- model$survey.design$variables -if (factor1_name %in% names(model_data_vars) && is.factor(model_data_vars[[factor1_name]])) { -levels1 <- levels(model_data_vars[[factor1_name]]) -} -if (factor2_name %in% names(model_data_vars) && is.factor(model_data_vars[[factor2_name]])) { -levels2 <- levels(model_data_vars[[factor2_name]]) -} -} else if (!is.null(model$data)) { -# For some other model types that might store data -if (factor1_name %in% names(model$data) && is.factor(model$data[[factor1_name]])) { -levels1 <- levels(model$data[[factor1_name]]) -} -if (factor2_name %in% names(model$data) && is.factor(model$data[[factor2_name]])) { -levels2 <- levels(model$data[[factor2_name]]) -} -} -if (is.null(levels1) || length(levels1) < 2) { -warning("Could not automatically determine levels for factor1: ", factor1_name, ". Check if it's a factor with >1 level in the model data.") -return(NULL) -} -if (is.null(levels2) || length(levels2) < 2) { -warning("Could not automatically determine levels for factor2: ", factor2_name, ". Check if it's a factor with >1 level in the model data.") -return(NULL) -} -ref_level1 <- levels1[1] -ref_level2 <- levels2[1] -non_ref_levels1 <- levels1[-1] -non_ref_levels2 <- levels2[-1] -all_results <- list() -# --- Loop through combinations of non-reference levels --- -for (level1 in non_ref_levels1) { -for (level2 in non_ref_levels2) { -comparison_name <- paste0(factor1_name, "_", level1, "_vs_", factor2_name, "_", level2) -message("Calculating for: ", comparison_name) -# --- Construct coefficient names --- -exp1_coef_name <- paste0(factor1_name, level1) -exp2_coef_name <- paste0(factor2_name, level2) -inter_coef_name1 <- paste0(factor1_name, level1, ":", factor2_name, level2) -inter_coef_name2 <- paste0(factor2_name, level2, ":", factor1_name, level1) -inter_coef_name <- NULL -if (inter_coef_name1 %in% coef_names_all) { -inter_coef_name <- inter_coef_name1 -} else if (inter_coef_name2 %in% coef_names_all) { -inter_coef_name <- inter_coef_name2 -} -# Check if all needed coefficients were found -current_coef_names_list <- list( -exp1_coef = exp1_coef_name, -exp2_coef = exp2_coef_name, -inter_coef = inter_coef_name -) -required_in_model <- c(exp1_coef_name, exp2_coef_name, inter_coef_name) -if (is.null(inter_coef_name) || any(!required_in_model %in% coef_names_all) ) { -warning("Could not find all necessary coefficients for interaction between ", -level1, " and ", level2, ". Skipping this combination.") -all_results[[comparison_name]] <- list(error = "Coefficients not found") -next # Skip to the next combination -} -# --- Call the calculation function 'addint' --- -interaction_results <- addint( # Call the dependency function -model = model, -type = "interaction", -coef_names = current_coef_names_list, -measures = measures, -conf.level = conf.level -) -all_results[[comparison_name]] <- interaction_results -} -} -return(all_results) -} -# Ensure the 'addint' function is defined or loaded from the package -# source("R/addint.R") # If running interactively before package build -#' -if (requireNamespace("survey", quietly = TRUE) && -requireNamespace("NHANES", quietly = TRUE) && -requireNamespace("dplyr", quietly = TRUE) && -requireNamespace("tidyr", quietly = TRUE) && -requireNamespace("msm", quietly = TRUE)) { -#' -library(survey) -library(NHANES) -library(dplyr) -library(tidyr) -library(msm) -#' -# --- 1. Data Preparation (NHANES Example) --- -data(NHANESraw) -#' -vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -"SDMVPSU", "SDMVSTRA", "WTMEC2YR") -#' -nhanes_adults_processed <- NHANESraw %>% -filter(Age >= 20) %>% -mutate( -ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -levels = c("Not Obese", "Obese")), -Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -levels = c("No", "Yes")), -Race1 = relevel(as.factor(Race1), ref = "White") -) %>% -select(all_of(vars_needed)) %>% -drop_na() -#' -adult_design_binary <- svydesign( -id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -nest = TRUE, data = nhanes_adults_processed -) -#' -# --- 2. Fit Interaction Term Model --- -interaction_model_logit <- svyglm( -Hypertension_130 ~ Race1 * ObeseStatus + Age, -design = adult_design_binary, family = quasibinomial() -) -#' -# --- 3. Calculate all additive interactions --- -all_interactions_list <- addintlist( -model = interaction_model_logit, -factor1_name = "Race1", -factor2_name = "ObeseStatus", -measures = "all" -) -#' -# Print the results list -print(all_interactions_list) -#' -# Example: Extract RERI for a specific interaction -# print(all_interactions_list$Race1_Black_vs_ObeseStatus_Obese$RERI) -#' -} else { -print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -} -} -#' Calculate All Additive Interaction Measures from an Interaction Model -#' -#' @description -#' Automatically identifies interaction terms in a fitted model involving two -#' categorical factors (e.g., `~ A * B`) and calculates measures of additive -#' interaction (RERI, AP, S) for each combination of non-reference levels. -#' This function requires the `addint()` function to be available. -#' -#' @details -#' This function serves as a wrapper around `addint()`. It determines the -#' levels of the two specified factors from the model object, constructs the -#' expected coefficient names for main effects and interaction terms (checking -#' both `A:B` and `B:A` formats), and then calls `addint()` with -#' `type = "interaction"` for each valid combination found in the model. -#' The results are compiled into a single summary table. -#' -#' @param model A fitted model object (e.g., `svycoxph`, `svyglm`) of type -#' `"interaction"` (containing a `factor1 * factor2` term). -#' @param factor1_name Character string: The name of the first factor variable -#' in the interaction (e.g., `"Race1"`). -#' @param factor2_name Character string: The name of the second factor variable -#' (e.g., `"ObeseStatus"`). -#' @param measures A character vector specifying which measures to calculate via `addint()`. -#' Options: `"RERI"`, `"AP"`, `"S"`, or `"all"`. Default is `"all"`. -#' @param conf.level Confidence level for the interval (default 0.95). -#' @param digits Integer. Number of decimal places for rounding estimates and -#' CIs in the final table. Default is 3. (Set to NULL for no rounding). -#' -#' @return A `tibble` (data frame) summarizing the additive interaction results. -#' Columns include: -#' \itemize{ -#' \item `Factor1`: Name of the first factor. -#' \item `Level1`: Non-reference level of the first factor for the comparison. -#' \item `Factor2`: Name of the second factor. -#' \item `Level2`: Non-reference level of the second factor for the comparison. -#' \item `Measure`: The additive interaction measure calculated ("RERI", "AP", or "S"). -#' \item `Estimate`: Point estimate of the measure. -#' \item `SE`: Standard Error (Note: For S, this is SE of log(S)). -#' \item `CI_low`: Lower confidence interval bound. -#' \item `CI_upp`: Upper confidence interval bound. -#' } -#' Returns an empty tibble if errors occur. -#' -#' @seealso \code{\link{addint}} -#' -#' @importFrom stats coef -#' @importFrom tibble tibble -#' @importFrom dplyr bind_rows mutate relocate select -#' @importFrom tidyr pivot_longer -#' @importFrom rlang := -#' -#' @export -#' -#' @examples -#' \donttest{ -#' # --- Load required libraries for the example --- -#' # Ensure the 'addint' function is defined or loaded from the package -#' # source("R/addint.R") # If running interactively before package build -#' -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE) && -#' requireNamespace("tidyr", quietly = TRUE) && -#' requireNamespace("msm", quietly = TRUE)) { -#' -#' library(survey) -#' library(NHANES) -#' library(dplyr) -#' library(tidyr) -#' library(msm) -#' -#' # --- 1. Data Preparation (NHANES Example) --- -#' data(NHANESraw) -#' -#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR") -#' -#' nhanes_adults_processed <- NHANESraw %>% -#' filter(Age >= 20) %>% -#' mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")), -#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -#' levels = c("No", "Yes")), -#' Race1 = relevel(as.factor(Race1), ref = "White") -#' ) %>% -#' select(all_of(vars_needed)) %>% -#' drop_na() -#' -#' adult_design_binary <- svydesign( -#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -#' nest = TRUE, data = nhanes_adults_processed -#' ) -#' -#' # --- 2. Fit Interaction Term Model --- -#' interaction_model_logit <- svyglm( -#' Hypertension_130 ~ Race1 * ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 3. Calculate all additive interactions --- -#' all_interactions_table <- addintlist( -#' model = interaction_model_logit, -#' factor1_name = "Race1", -#' factor2_name = "ObeseStatus", -#' measures = "all" -#' ) -#' -#' # Print the results table -#' print(all_interactions_table, n=50) -#' -#' } else { -#' print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -#' } -#' } -addintlist <- function(model, -factor1_name, -factor2_name, -measures = "all", -conf.level = 0.95, -digits = 3) { # Added digits argument -# Check if the core calculation function 'addint' exists -if (!exists("addint") || !is.function(addint)) { -stop("The 'addint' function is required but was not found. Ensure it is loaded.", call. = FALSE) -} -valid_measures_req <- c("RERI", "AP", "S", "all") -if (any(!measures %in% valid_measures_req)) { -stop("Invalid measure specified in 'measures'. Choose from 'RERI', 'AP', 'S', or 'all'.", call. = FALSE) -} -if ("all" %in% measures) { -measures_to_calc <- c("RERI", "AP", "S") -} else { -measures_to_calc <- unique(measures) -} -if (!inherits(model, c("svycoxph", "svyglm", "coxph", "glm"))) { -warning("Model object may not be of a supported type.") -# Attempt to proceed anyway -} -beta <- tryCatch(stats::coef(model), error = function(e) NULL) -if (is.null(beta)) { -warning("Could not extract coefficients from the model.") -return(tibble::tibble()) # Return empty tibble -} -coef_names_all <- names(beta) -# --- Get Factor Levels (excluding reference) --- -levels1 <- NULL -levels2 <- NULL -# Try accessing xlevels if they exist (common in lm, glm) -if (!is.null(model$xlevels)) { -levels1 <- model$xlevels[[factor1_name]] -levels2 <- model$xlevels[[factor2_name]] -} else if (inherits(model, "survey.design")) { -# For survey design objects passed directly (less common) -if (factor1_name %in% names(model$variables) && is.factor(model$variables[[factor1_name]])) { -levels1 <- levels(model$variables[[factor1_name]]) -} -if (factor2_name %in% names(model$variables) && is.factor(model$variables[[factor2_name]])) { -levels2 <- levels(model$variables[[factor2_name]]) -} -} else if (!is.null(model$survey.design)) { -# For survey model objects (svyglm, svycoxph) -model_data_vars <- model$survey.design$variables -if (factor1_name %in% names(model_data_vars) && is.factor(model_data_vars[[factor1_name]])) { -levels1 <- levels(model_data_vars[[factor1_name]]) -} -if (factor2_name %in% names(model_data_vars) && is.factor(model_data_vars[[factor2_name]])) { -levels2 <- levels(model_data_vars[[factor2_name]]) -} -} else if (!is.null(model$data)) { -# For some other model types that might store data -if (factor1_name %in% names(model$data) && is.factor(model$data[[factor1_name]])) { -levels1 <- levels(model$data[[factor1_name]]) -} -if (factor2_name %in% names(model$data) && is.factor(model$data[[factor2_name]])) { -levels2 <- levels(model$data[[factor2_name]]) -} -} -if (is.null(levels1) || length(levels1) < 2) { -warning("Could not automatically determine levels for factor1: ", factor1_name, ". Check if it's a factor with >1 level in the model data.") -return(tibble::tibble()) # Return empty tibble -} -if (is.null(levels2) || length(levels2) < 2) { -warning("Could not automatically determine levels for factor2: ", factor2_name, ". Check if it's a factor with >1 level in the model data.") -return(tibble::tibble()) # Return empty tibble -} -ref_level1 <- levels1[1] -ref_level2 <- levels2[1] -non_ref_levels1 <- levels1[-1] -non_ref_levels2 <- levels2[-1] -all_results_list <- list() # Store results before binding -# --- Loop through combinations of non-reference levels --- -for (level1 in non_ref_levels1) { -for (level2 in non_ref_levels2) { -comparison_name <- paste0(factor1_name, "_", level1, "_vs_", factor2_name, "_", level2) -message("Calculating for: ", comparison_name) -# --- Construct coefficient names --- -exp1_coef_name <- paste0(factor1_name, level1) -exp2_coef_name <- paste0(factor2_name, level2) -inter_coef_name1 <- paste0(factor1_name, level1, ":", factor2_name, level2) -inter_coef_name2 <- paste0(factor2_name, level2, ":", factor1_name, level1) -inter_coef_name <- NULL -if (inter_coef_name1 %in% coef_names_all) { -inter_coef_name <- inter_coef_name1 -} else if (inter_coef_name2 %in% coef_names_all) { -inter_coef_name <- inter_coef_name2 -} -# Check if all needed coefficients were found -current_coef_names_list <- list( -exp1_coef = exp1_coef_name, -exp2_coef = exp2_coef_name, -inter_coef = inter_coef_name -) -required_in_model <- c(exp1_coef_name, exp2_coef_name, inter_coef_name) -if (is.null(inter_coef_name) || any(!required_in_model %in% coef_names_all) ) { -warning("Could not find all necessary coefficients for interaction between ", -level1, " and ", level2, ". Skipping this combination.") -# Store minimal info about the skip -all_results_list[[length(all_results_list) + 1]] <- tibble::tibble( -Factor1 = factor1_name, Level1 = level1, -Factor2 = factor2_name, Level2 = level2, -Measure = "Error", Estimate = NA_real_, SE = NA_real_, -CI_low = NA_real_, CI_upp = NA_real_ -) -next # Skip to the next combination -} -# --- Call the calculation function 'addint' --- -interaction_results <- addint( # Call the dependency function -model = model, -type = "interaction", -coef_names = current_coef_names_list, -measures = measures_to_calc, # Use the validated list -conf.level = conf.level -) -# --- Process results into a flat tibble for this combination --- -if (!is.null(interaction_results) && length(interaction_results) > 0) { -for(measure_name in names(interaction_results)) { -res_vector <- interaction_results[[measure_name]] -se_val <- ifelse(measure_name == "S", res_vector["S_SE_log"], res_vector[2]) # Get correct SE column -all_results_list[[length(all_results_list) + 1]] <- tibble::tibble( -Factor1 = factor1_name, -Level1 = level1, -Factor2 = factor2_name, -Level2 = level2, -Measure = measure_name, -Estimate = res_vector[1], # First element is always estimate -SE = se_val, -CI_low = res_vector[3], # Third element is always LowerCI -CI_upp = res_vector[4] # Fourth element is always UpperCI -) -} -} else { -# Store minimal info if addint returned NULL -all_results_list[[length(all_results_list) + 1]] <- tibble::tibble( -Factor1 = factor1_name, Level1 = level1, -Factor2 = factor2_name, Level2 = level2, -Measure = "Error", Estimate = NA_real_, SE = NA_real_, -CI_low = NA_real_, CI_upp = NA_real_ -) -warning("Calculation failed within addint() for interaction between ", -level1, " and ", level2, ". Check previous warnings.") -} -} # end loop level2 -} # end loop level1 -# --- Combine all results into a single tibble --- -if (length(all_results_list) == 0) { -warning("No interaction results could be calculated.") -# Return an empty tibble with correct columns -return(tibble::tibble(Factor1 = character(), Level1 = character(), -Factor2 = character(), Level2 = character(), -Measure = character(), Estimate = double(), SE = double(), -CI_low = double(), CI_upp = double())) -} -final_table <- dplyr::bind_rows(all_results_list) -# --- Optional rounding --- -if (!is.null(digits)) { -final_table <- final_table %>% -dplyr::mutate(dplyr::across(c(Estimate, SE, CI_low, CI_upp), ~round(.x, digits))) -} -# --- Final structure --- -final_table <- final_table %>% -dplyr::relocate(Factor1, Level1, Factor2, Level2, Measure, Estimate, SE, CI_low, CI_upp) -return(final_table) -} -# --- Load required libraries for the example --- -# Ensure the 'addint' function is defined or loaded from the package -# source("R/addint.R") # If running interactively before package build -#' -if (requireNamespace("survey", quietly = TRUE) && -requireNamespace("NHANES", quietly = TRUE) && -requireNamespace("dplyr", quietly = TRUE) && -requireNamespace("tidyr", quietly = TRUE) && -requireNamespace("msm", quietly = TRUE)) { -#' -library(survey) -library(NHANES) -library(dplyr) -library(tidyr) -library(msm) -#' -# --- 1. Data Preparation (NHANES Example) --- -data(NHANESraw) -#' -vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -"SDMVPSU", "SDMVSTRA", "WTMEC2YR") -#' -nhanes_adults_processed <- NHANESraw %>% -filter(Age >= 20) %>% -mutate( -ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -levels = c("Not Obese", "Obese")), -Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -levels = c("No", "Yes")), -Race1 = relevel(as.factor(Race1), ref = "White") -) %>% -select(all_of(vars_needed)) %>% -drop_na() -#' -adult_design_binary <- svydesign( -id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -nest = TRUE, data = nhanes_adults_processed -) -#' -# --- 2. Fit Interaction Term Model --- -interaction_model_logit <- svyglm( -Hypertension_130 ~ Race1 * ObeseStatus + Age, -design = adult_design_binary, family = quasibinomial() -) -#' -# --- 3. Calculate all additive interactions --- -all_interactions_table <- addintlist( -model = interaction_model_logit, -factor1_name = "Race1", -factor2_name = "ObeseStatus", -measures = "all" -) -#' -# Print the results table -print(all_interactions_table, n=50) -#' -} else { -print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -} -} diff --git a/R/addint.R b/R/addint.R index b24f9af..255e0f7 100644 --- a/R/addint.R +++ b/R/addint.R @@ -48,98 +48,33 @@ #' #' @examples #' \donttest{ -#' # --- Load required libraries for the example --- -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE) && -#' requireNamespace("tidyr", quietly = TRUE) && -#' requireNamespace("msm", quietly = TRUE)) { +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' library(survey) -#' library(NHANES) -#' library(dplyr) -#' library(tidyr) -#' library(msm) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # --- 1. Data Preparation (NHANES Example) --- -#' data(NHANESraw) +#' # Interaction-term model: additive interaction of sex and insulin use. +#' fit <- survey::svyglm(htn01 ~ sex * insulin + age, +#' design = design, family = quasibinomial()) #' -#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR") +#' # Coefficient names must match names(coef(fit)) exactly. +#' addint( +#' model = fit, +#' type = "interaction", +#' coef_names = list( +#' exp1_coef = "sexFemale", +#' exp2_coef = "insulinYes", +#' inter_coef = "sexFemale:insulinYes" +#' ), +#' measures = "all" +#' ) #' -#' nhanes_adults_processed <- NHANESraw %>% -#' filter(Age >= 20) %>% -#' mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")), -#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -#' levels = c("No", "Yes")), -#' Race1 = relevel(as.factor(Race1), ref = "White") -#' ) %>% -#' select(all_of(vars_needed)) %>% -#' drop_na() -#' -#' adult_design_binary <- svydesign( -#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -#' nest = TRUE, data = nhanes_adults_processed -#' ) -#' -#' # --- 2a. Fit Joint Variable Model --- -#' adult_design_binary <- update(adult_design_binary, -#' Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE) -#' ) -#' adult_design_binary <- update(adult_design_binary, -#' Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese") -#' ) -#' -#' joint_model_logit <- svyglm( -#' Hypertension_130 ~ Race1_ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 2b. Fit Interaction Term Model --- -#' interaction_model_logit <- svyglm( -#' Hypertension_130 ~ Race1 * ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 3. Calculate Additive Interaction (Black vs White * Obese vs Not Obese) --- -#' -#' # --- Using the Joint Model --- -#' joint_coef_names_black_obese <- list( -#' exp1_level = "Race1_ObeseStatusBlack_Not Obese", # A=1, B=0 -#' exp2_level = "Race1_ObeseStatusWhite_Obese", # A=0, B=1 -#' both_levels = "Race1_ObeseStatusBlack_Obese" # A=1, B=1 -#' ) -#' -#' results_joint <- addint( -#' model = joint_model_logit, -#' type = "joint", -#' coef_names = joint_coef_names_black_obese, -#' measures = "all" -#' ) -#' print("--- Results from Joint Model (Black vs Obese Interaction) ---") -#' print(results_joint) -#' -#' # --- Using the Interaction Model --- -#' interaction_coef_names_black_obese <- list( -#' exp1_coef = "Race1Black", # Main effect A -#' exp2_coef = "ObeseStatusObese", # Main effect B -#' inter_coef = "Race1Black:ObeseStatusObese" # Interaction A:B -#' ) -#' -#' results_interact <- addint( -#' model = interaction_model_logit, -#' type = "interaction", -#' coef_names = interaction_coef_names_black_obese, -#' measures = "all" -#' ) -#' print("--- Results from Interaction Model (Black vs Obese Interaction) ---") -#' print(results_interact) -#' -#' } else { -#' print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -#' } +#' # RERI/AP/S can also be computed from a joint-variable model with +#' # type = "joint" (see the 'coef_names' argument). For logistic models these +#' # use odds ratios, which approximate risk ratios only when the outcome is rare. #' } addint <- function(model, type = c("joint", "interaction"), diff --git a/R/addintlist.R b/R/addintlist.R index 4964e7c..80c2743 100644 --- a/R/addintlist.R +++ b/R/addintlist.R @@ -25,6 +25,8 @@ #' @param conf.level Confidence level for the interval (default 0.95). #' @param digits Integer. Number of decimal places for rounding estimates and #' CIs in the final table. Default is 3. (Set to NULL for no rounding). +#' @param verbose Logical. If `TRUE`, prints a progress message for each +#' factor-level combination as it is calculated. Defaults to `FALSE`. #' #' @return A `tibble` (data frame) summarizing the additive interaction results. #' Columns include: @@ -118,7 +120,8 @@ addintlist <- function(model, factor2_name, measures = "all", conf.level = 0.95, - digits = 3) { # Added digits argument + digits = 3, + verbose = FALSE) { # Check if the core calculation function 'addint' exists if (!exists("addint") || !is.function(addint)) { @@ -205,7 +208,7 @@ addintlist <- function(model, for (level2 in non_ref_levels2) { comparison_name <- paste0(factor1_name, "_", level1, "_vs_", factor2_name, "_", level2) - message("Calculating for: ", comparison_name) + if (verbose) message("Calculating for: ", comparison_name) # --- Construct coefficient names --- exp1_coef_name <- paste0(factor1_name, level1) @@ -306,8 +309,8 @@ addintlist <- function(model, # --- Final structure --- final_table <- final_table %>% # Use .data pronoun - dplyr::relocate(.data$Factor1, .data$Level1, .data$Factor2, .data$Level2, - .data$Measure, .data$Estimate, .data$SE, .data$CI_low, .data$CI_upp) + dplyr::relocate("Factor1", "Level1", "Factor2", "Level2", + "Measure", "Estimate", "SE", "CI_low", "CI_upp") return(final_table) } diff --git a/R/inteffects.R b/R/inteffects.R index 9c4238f..c06b158 100644 --- a/R/inteffects.R +++ b/R/inteffects.R @@ -50,87 +50,33 @@ #' #' @examples #' \donttest{ -#' # --- Load required libraries for the example --- -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE) && -#' requireNamespace("tidyr", quietly = TRUE) && -#' requireNamespace("msm", quietly = TRUE)) { # msm needed for SE calculation +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' library(survey) -#' library(NHANES) -#' library(dplyr) -#' library(tidyr) -#' library(msm) # Load msm +#' # Build a joint variable from two factors, then fit a joint-variable model. +#' nhanes_mortality$sex_insulin <- interaction( +#' nhanes_mortality$sex, nhanes_mortality$insulin, sep = "_" +#' ) #' -#' # --- 1. Data Preparation (NHANES Example) --- -#' data(NHANESraw) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR") +#' joint_fit <- survey::svyglm(htn01 ~ sex_insulin + age, +#' design = design, family = quasibinomial()) #' -#' nhanes_adults_processed <- NHANESraw %>% -#' filter(Age >= 20) %>% -#' mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")), -#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -#' levels = c("No", "Yes")), -#' Race1 = relevel(as.factor(Race1), ref = "White") -#' ) %>% -#' select(all_of(vars_needed)) %>% -#' drop_na() -#' -#' adult_design_binary <- svydesign( -#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -#' nest = TRUE, data = nhanes_adults_processed -#' ) -#' -#' # --- 2. Create Joint Variable and Fit Joint Model --- -#' adult_design_binary <- update(adult_design_binary, -#' Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE) -#' ) -#' adult_design_binary <- update(adult_design_binary, -#' Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese") -#' ) -#' -#' joint_model_logit <- svyglm( -#' Hypertension_130 ~ Race1_ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 3. Run the function --- -#' f1_levels <- levels(adult_design_binary$variables$Race1) -#' f2_levels <- levels(adult_design_binary$variables$ObeseStatus) -#' -#' # --- Example 1: Output on Ratio scale --- -#' simple_effects_ratio <- inteffects( -#' joint_model = joint_model_logit, -#' joint_var_name = "Race1_ObeseStatus", -#' factor1_name = "Race1", -#' factor2_name = "ObeseStatus", -#' factor1_levels = f1_levels, -#' factor2_levels = f2_levels, -#' level_separator = "_", -#' scale = "ratio" -#' ) -#' print("--- Output on Ratio Scale ---") -#' print(simple_effects_ratio, n = 50) -#' -#' # --- Example 2: Output on Log scale --- -#' simple_effects_log <- inteffects( -#' joint_model = joint_model_logit, -#' joint_var_name = "Race1_ObeseStatus", -#' factor1_name = "Race1", -#' factor2_name = "ObeseStatus", -#' factor1_levels = f1_levels, -#' factor2_levels = f2_levels, -#' level_separator = "_", -#' scale = "log" -#' ) -#' print("--- Output on Log Scale ---") -#' print(simple_effects_log, n = 50) -#' } +#' # Stratum-specific (simple) effects reconstructed from the joint model. +#' inteffects( +#' joint_model = joint_fit, +#' joint_var_name = "sex_insulin", +#' factor1_name = "sex", +#' factor2_name = "insulin", +#' factor1_levels = c("Male", "Female"), # reference level first +#' factor2_levels = c("No", "Yes"), # reference level first +#' level_separator = "_", +#' scale = "ratio" +#' ) #' } inteffects <- function(joint_model, joint_var_name, @@ -322,12 +268,12 @@ inteffects <- function(joint_model, output_df <- output_df %>% # Use .data pronoun dplyr::select( - .data$Comparison, - Estimate = .data$logEstimate, - SE = .data$SE_log, - CI.low = .data$CI.low_log, - CI.upp = .data$CI.upp_log, - .data$p_value + "Comparison", + Estimate = "logEstimate", + SE = "SE_log", + CI.low = "CI.low_log", + CI.upp = "CI.upp_log", + "p_value" ) } else { # scale == "ratio" @@ -346,12 +292,12 @@ inteffects <- function(joint_model, output_df <- output_df %>% # Use .data pronoun dplyr::select( - .data$Comparison, - Estimate = .data$Estimate_ratio, - SE = .data$SE_ratio, # Use ratio SE - CI.low = .data$CI.low_ratio, - CI.upp = .data$CI.upp_ratio, - .data$p_value + "Comparison", + Estimate = "Estimate_ratio", + SE = "SE_ratio", # Use ratio SE + CI.low = "CI.low_ratio", + CI.upp = "CI.upp_ratio", + "p_value" ) } @@ -359,10 +305,9 @@ inteffects <- function(joint_model, output_df <- output_df %>% # Use .data pronoun dplyr::rename( - `p-value` = .data$p_value + `p-value` = "p_value" ) %>% - # Use .data pronoun - dplyr::relocate(.data$Comparison, .data$Estimate, .data$SE, .data$CI.low, .data$CI.upp, .data$`p-value`) + dplyr::relocate("Comparison", "Estimate", "SE", "CI.low", "CI.upp", "p-value") return(output_df) diff --git a/R/jointeffects.R b/R/jointeffects.R index 6cc3958..3060f9e 100644 --- a/R/jointeffects.R +++ b/R/jointeffects.R @@ -346,12 +346,12 @@ jointeffects <- function(interaction_model, final_df <- output_df %>% # Use .data pronoun dplyr::select( - .data$Level1, - .data$Level2, - Estimate = .data$logEstimate, - SE = .data$SE_log, - CI.low = .data$CI.low_log, - CI.upp = .data$CI.upp_log + "Level1", + "Level2", + Estimate = "logEstimate", + SE = "SE_log", + CI.low = "CI.low_log", + CI.upp = "CI.upp_log" ) } else { # scale == "ratio" # Round ratio scale columns if digits specified @@ -366,12 +366,12 @@ jointeffects <- function(interaction_model, final_df <- output_df %>% # Use .data pronoun dplyr::select( - .data$Level1, - .data$Level2, - Estimate = .data$Estimate_ratio, - SE = .data$SE_ratio, - CI.low = .data$CI.low_ratio, - CI.upp = .data$CI.upp_ratio + "Level1", + "Level2", + Estimate = "Estimate_ratio", + SE = "SE_ratio", + CI.low = "CI.low_ratio", + CI.upp = "CI.upp_ratio" # CI_formatted column is now completely removed ) } diff --git a/R/reportint.R b/R/reportint.R index f9d2934..d25b91b 100644 --- a/R/reportint.R +++ b/R/reportint.R @@ -57,71 +57,28 @@ #' @export #' #' @examples -#' \dontrun{ -#' # --- Load required libraries for the example --- -#' if (requireNamespace("svyTable1", quietly = TRUE) && -#' requireNamespace("Publish", quietly = TRUE) && -#' requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE) && -#' requireNamespace("tidyr", quietly = TRUE) && -#' requireNamespace("broom", quietly = TRUE)) { +#' \donttest{ +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' library(svyTable1) -#' library(Publish) -#' library(survey) -#' library(NHANES) -#' library(dplyr) -#' library(tidyr) -#' library(broom) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # --- 1. Data Preparation (NHANES Example) --- -#' data(NHANESraw) -#' -#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", -#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR") -#' -#' nhanes_adults_processed <- NHANESraw %>% -#' dplyr::filter(Age >= 20) %>% -#' dplyr::mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")), -#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), -#' levels = c("No", "Yes")), -#' Race1 = relevel(as.factor(Race1), ref = "White") -#' ) %>% -#' dplyr::select(dplyr::all_of(vars_needed)) %>% -#' tidyr::drop_na() -#' -#' adult_design_binary <- svydesign( -#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -#' nest = TRUE, data = nhanes_adults_processed -#' ) -#' -#' # --- 2. FIT ONLY THE INTERACTION MODEL --- -#' -#' interaction_model_fit <- svyglm( -#' Hypertension_130 ~ Race1 * ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 3. Run the wrapper function --- -#' -#' report_list <- reportint( -#' interaction_model = interaction_model_fit -#' ) -#' -#' # You can then print the individual panel tables: -#' # print(report_list$joint_effects) -#' # print(report_list$stratum_specific_effects) -#' -#' # --- OR --- -#' -#' # Print the single, publication-ready summary tables: -#' # print(report_list$effect_modification_report) -#' # print(report_list$Interaction_report) +#' # Fit a model with a two-factor interaction (sex x insulin). +#' interaction_fit <- survey::svyglm( +#' htn01 ~ sex * insulin + age, +#' design = design, family = quasibinomial() +#' ) #' -#' } +#' # The default output = "list" returns the panel data frames (no rendering). +#' report <- reportint(interaction_fit, +#' factor1_name = "sex", factor2_name = "insulin", +#' output = "list") +#' report$joint_effects +#' report$additive_interaction +#' report$Interaction_report #' } reportint <- function(interaction_model, factor1_name = NULL, @@ -163,13 +120,13 @@ reportint <- function(interaction_model, # --- For "viewer" or "rmd", we must render the Rmd template --- - # Find the Rmd template within your package - # *** YOU MUST REPLACE "YourPackageName" WITH YOUR ACTUAL PACKAGE NAME *** + # Find the Rmd template that ships with the installed package. template_path <- system.file("rmd/interaction_report.Rmd", package = "svyTable1") if (template_path == "") { - stop("Could not find 'interaction_report.Rmd' in 'inst/rmd/'.\n", - "Make sure the file exists and 'YourPackageName' is correct.", call. = FALSE) + stop("Could not find the interaction report template ", + "(inst/rmd/interaction_report.Rmd) in the installed 'svyTable1' package. ", + "Try reinstalling the package, or use output = \"list\".", call. = FALSE) } # --- Create render parameters --- @@ -348,7 +305,7 @@ reportint <- function(interaction_model, ) panelA_table <- panelA_data_ratio %>% - dplyr::left_join(dplyr::select(panelA_data_log, .data$Level1, .data$Level2, .data$p_value), by = c("Level1", "Level2")) %>% + dplyr::left_join(dplyr::select(panelA_data_log, "Level1", "Level2", "p_value"), by = c("Level1", "Level2")) %>% dplyr::mutate( OR_str = sprintf(paste0("%.", digits_ratio, "f"), .data$Estimate), CI_str = format_ci(.data$CI.low, .data$CI.upp, digits = digits_ratio), @@ -357,11 +314,11 @@ reportint <- function(interaction_model, CI_str = dplyr::if_else(.data$Level1 == ref_level1 & .data$Level2 == ref_level2, "(Reference)", .data$CI_str) ) %>% dplyr::select( - !!rlang::sym(factor1_name) := .data$Level1, - !!rlang::sym(factor2_name) := .data$Level2, - OR = .data$OR_str, - `95% CI` = .data$CI_str, - p = .data$p_str + !!rlang::sym(factor1_name) := "Level1", + !!rlang::sym(factor2_name) := "Level2", + OR = "OR_str", + `95% CI` = "CI_str", + p = "p_str" ) # --- 6. Generate & Format Panel B --- @@ -369,7 +326,11 @@ reportint <- function(interaction_model, ci_format_str <- "(%s, %s)" - pub_obj <- Publish::publish(interaction_model, ci.format = ci_format_str, digits = digits_ratio) + # Publish::publish() emits a deprecation note about its internal 'robust' + # argument (from emmeans); it is unrelated to our call, so silence it. + pub_obj <- suppressWarnings( + Publish::publish(interaction_model, ci.format = ci_format_str, digits = digits_ratio) + ) panelB_table_raw <- pub_obj$regressionTable table_names <- names(panelB_table_raw) @@ -390,12 +351,12 @@ reportint <- function(interaction_model, panelB_table <- panelB_table_raw %>% dplyr::filter(grepl(factor1_name, .data$Variable, ignore.case = TRUE) | grepl(factor2_name, .data$Variable, ignore.case = TRUE)) %>% dplyr::rename( - Comparison = .data$Variable, + Comparison = "Variable", Estimate = !!rlang::sym(estimate_col_name), - `95% CI` = .data$CI.95, - p = .data$`p-value` + `95% CI` = "CI.95", + p = "p-value" ) %>% - dplyr::select(.data$Comparison, .data$Estimate, .data$`95% CI`, .data$p) + dplyr::select("Comparison", "Estimate", "95% CI", "p") # --- 7. Generate & Format Panel C --- @@ -406,7 +367,8 @@ reportint <- function(interaction_model, factor2_name = factor2_name, measures = "all", conf.level = conf.level, - digits = digits_additive + digits = digits_additive, + verbose = verbose ) panelC_table <- panelC_data_raw %>% @@ -415,21 +377,21 @@ reportint <- function(interaction_model, Est_str = sprintf(paste0("%.", digits_additive, "f"), .data$Estimate), CI_str = format_ci(.data$CI_low, .data$CI_upp, digits = digits_additive) ) %>% - dplyr::select(.data$Level1, .data$Level2, .data$Measure, .data$Est_str, .data$CI_str) %>% + dplyr::select("Level1", "Level2", "Measure", "Est_str", "CI_str") %>% tidyr::pivot_wider( - names_from = .data$Measure, - values_from = c(.data$Est_str, .data$CI_str), + names_from = "Measure", + values_from = c("Est_str", "CI_str"), names_glue = "{Measure}_{.value}" ) %>% dplyr::select( - !!rlang::sym(factor1_name) := .data$Level1, - !!rlang::sym(factor2_name) := .data$Level2, - RERI = .data$RERI_Est_str, - `RERI 95% CI` = .data$RERI_CI_str, - AP = .data$AP_Est_str, - `AP 95% CI` = .data$AP_CI_str, - S = .data$S_Est_str, - `S 95% CI` = .data$S_CI_str + !!rlang::sym(factor1_name) := "Level1", + !!rlang::sym(factor2_name) := "Level2", + RERI = "RERI_Est_str", + `RERI 95% CI` = "RERI_CI_str", + AP = "AP_Est_str", + `AP 95% CI` = "AP_CI_str", + S = "S_Est_str", + `S 95% CI` = "S_CI_str" ) # --- 8. Generate Multiplicative Interaction --- @@ -457,10 +419,10 @@ reportint <- function(interaction_model, p_str = format_p_value(.data$p.value) ) %>% dplyr::select( - InteractionTerm = .data$term, - ROR = .data$ROR_str, - `95% CI` = .data$CI_str, - p = .data$p_str + InteractionTerm = "term", + ROR = "ROR_str", + `95% CI` = "CI_str", + p = "p_str" ) @@ -582,27 +544,27 @@ reportint <- function(interaction_model, ) %>% dplyr::select( "Characteristic", - Estimate = .data$OR, - `95% CI` = .data$`95% CI`, - `p-value` = .data$p + Estimate = "OR", + `95% CI` = "95% CI", + `p-value` = "p" ) # 2. Format Panel B panelB_fmt <- panelB %>% dplyr::select( - Characteristic = .data$Comparison, - Estimate = .data$Estimate, - `95% CI` = .data$`95% CI`, - `p-value` = .data$p + Characteristic = "Comparison", + Estimate = "Estimate", + `95% CI` = "95% CI", + `p-value` = "p" ) # 3. Format Panel M (Multiplicative) panelM_fmt <- panelM %>% dplyr::select( - Characteristic = .data$InteractionTerm, - Estimate = .data$ROR, - `95% CI` = .data$`95% CI`, - `p-value` = .data$p + Characteristic = "InteractionTerm", + Estimate = "ROR", + `95% CI` = "95% CI", + `p-value` = "p" ) # 4. Format Panel C (Additive) @@ -612,12 +574,12 @@ reportint <- function(interaction_model, ) %>% dplyr::select( "Characteristic", - RERI = .data$RERI, - `RERI 95% CI` = .data$`RERI 95% CI`, - AP = .data$AP, - `AP 95% CI` = .data$`AP 95% CI`, - S = .data$S, - `S 95% CI` = .data$`S 95% CI` + RERI = "RERI", + `RERI 95% CI` = "RERI 95% CI", + AP = "AP", + `AP 95% CI` = "AP 95% CI", + S = "S", + `S 95% CI` = "S 95% CI" ) # 5. Create Header Rows @@ -649,7 +611,7 @@ reportint <- function(interaction_model, if(estimate_col_name == "Coefficient") est_col_name_final <- "Estimate" final_table <- final_table %>% - dplyr::rename(!!est_col_name_final := .data$Estimate) + dplyr::rename(!!est_col_name_final := "Estimate") return(final_table) } diff --git a/R/svyAUC.R b/R/svyAUC.R index 8800578..558444f 100644 --- a/R/svyAUC.R +++ b/R/svyAUC.R @@ -10,54 +10,38 @@ #' @param plot A logical value. If `TRUE`, an ROC curve is plotted. Defaults to `FALSE`. #' #' @return -#' If `plot = FALSE` (default), returns a `data.frame` with the AUC, SE, and 95% CI. -#' If `plot = TRUE`, invisibly returns a list containing the summary `data.frame` -#' and another `data.frame` with the ROC curve coordinates (TPR and FPR). +#' If `plot = FALSE` (default), a one-row `data.frame` with columns `AUC`, `SE`, +#' `CI_Lower` and `CI_Upper`. The standard error is computed across the design's +#' replicate weights and the 95\% confidence interval is built on the logit scale +#' (so it stays within `[0, 1]`); both are therefore conditional on the +#' replication scheme used to build `design`. If `plot = TRUE`, an ROC curve is +#' drawn and a list is returned invisibly with the summary `data.frame` and a +#' `data.frame` of ROC coordinates (TPR and FPR). #' -#' @importFrom survey withReplicates SE +#' @importFrom survey withReplicates SE degf #' @importFrom WeightedROC WeightedROC WeightedAUC -#' @importFrom stats model.frame model.matrix coef plogis weights +#' @importFrom stats model.frame model.matrix coef plogis qlogis qt weights #' @importFrom graphics plot abline title #' #' @export #' #' @examples -#' \dontrun{ -#' # Ensure required packages are loaded -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE)) { +#' \donttest{ +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' # 1. Prepare Data -#' data(NHANESraw, package = "NHANES") -#' nhanes_data <- NHANESraw %>% -#' dplyr::filter(Age >= 20) %>% -#' dplyr::mutate(ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese"))) %>% -#' dplyr::filter(complete.cases(ObeseStatus, Age, Gender, Race1, -#' WTMEC2YR, SDMVPSU, SDMVSTRA)) +#' # svyAUC requires a replicate-weights design. +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) +#' rep_design <- survey::as.svrepdesign(design) #' -#' # 2. Create a replicate design object -#' std_design <- survey::svydesign( -#' ids = ~SDMVPSU, -#' strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, -#' nest = TRUE, -#' data = nhanes_data -#' ) -#' rep_design <- survey::as.svrepdesign(std_design) +#' fit <- survey::svyglm(htn01 ~ age + sex + smoking, +#' design = rep_design, family = quasibinomial()) #' -#' # 3. Fit a survey logistic regression model using the replicate design -#' fit_obesity_rep <- survey::svyglm( -#' ObeseStatus ~ Age + Gender + Race1, -#' design = rep_design, -#' family = quasibinomial() -#' ) -#' -#' # 4. Calculate the design-correct AUC -#' auc_results <- svyAUC(fit_obesity_rep, rep_design) -#' print(auc_results) -#' } +#' # AUC with a design-correct SE and a logit-scale 95% CI. +#' svyAUC(fit, rep_design) #' } svyAUC <- function(fit, design, plot = FALSE) { @@ -108,25 +92,38 @@ svyAUC <- function(fit, design, plot = FALSE) { return.replicates = TRUE ) - # Manually calculate the confidence interval - auc_estimate <- result$theta - se <- survey::SE(result) + # Point estimate and design-correct SE (variance across replicate weights). + auc_estimate <- as.numeric(result$theta) + se <- as.numeric(survey::SE(result)) + + # Critical value from the design degrees of freedom rather than a fixed 1.96. + df <- survey::degf(design) + crit <- stats::qt(0.975, df) + + # Build the 95% CI on the logit scale and back-transform, so the interval + # cannot fall outside the valid [0, 1] range of an AUC. At the boundaries + # (where the logit is undefined) fall back to a clamped Wald interval. + if (is.finite(auc_estimate) && auc_estimate > 0 && auc_estimate < 1) { + se_logit <- se / (auc_estimate * (1 - auc_estimate)) + ci_lower <- stats::plogis(stats::qlogis(auc_estimate) - crit * se_logit) + ci_upper <- stats::plogis(stats::qlogis(auc_estimate) + crit * se_logit) + } else { + ci_lower <- max(0, auc_estimate - crit * se) + ci_upper <- min(1, auc_estimate + crit * se) + } summary_df <- data.frame( AUC = auc_estimate, SE = se, - CI_Lower = auc_estimate - 1.96 * se, - CI_Upper = auc_estimate + 1.96 * se + CI_Lower = ci_lower, + CI_Upper = ci_upper ) rownames(summary_df) <- NULL # --- PLOTTING LOGIC --- if (plot) { - # Calculate ROC curve points using the full-sample weights + # Calculate ROC curve points using the full-sample weights. full_weights <- weights(design, "sampling") - roc_data <- auc_statistic(full_weights, design$variables) # Temporarily re-run to get roc_curve - - # Actually need the curve, not just the AUC predictions <- predict(fit, newdata = design$variables, type = "response") outcome <- design$variables[[outcome_name]] if(is.factor(outcome)) { diff --git a/R/svycoxph_CE.R b/R/svycoxph_CE.R index a9f0033..b958feb 100644 --- a/R/svycoxph_CE.R +++ b/R/svycoxph_CE.R @@ -47,7 +47,10 @@ #' before splitting. If this fails, the user can manually provide the #' \code{design_ids}, \code{design_weights}, and \code{design_strata} arguments. #' -#' @return A \code{ggplot} object showing the coefficient trend over time. +#' @return A \code{ggplot} object. The y-axis is the estimated log hazard ratio +#' for \code{var_to_test} within each follow-up interval, plotted against time; +#' a roughly flat trend supports the proportional-hazards (constant-effect) +#' assumption. This is a visual diagnostic, not a formal hypothesis test. #' #' @importFrom survey svycoxph svydesign #' @importFrom survival survSplit @@ -59,52 +62,25 @@ #' @export #' #' @examples -#' \dontrun{ -#' # --- 1. Load data and create a design --- -#' library(survey) -#' library(dplyr) +#' \donttest{ #' data(nhanes_mortality, package = "svyTable1") #' -#' # Create a base design -#' analytic_design <- svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = nhanes_mortality, -#' nest = TRUE +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality #' ) #' -#' # Prepare data (filter to >0 time) -#' data_clean <- analytic_design$variables %>% -#' dplyr::filter(stime > 0) %>% -#' mutate(caff_bin = ifelse(caff == "No consumption", "No", "Yes")) -#' -#' # Create final design object -#' final_design <- svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = data_clean, -#' nest = TRUE -#' ) -#' -#' # --- Example 1: Standard Usage --- +#' # Visual check of the constant-effect (proportional hazards) assumption for +#' # 'age': a roughly flat trend across follow-up intervals supports it. #' svycoxph_CE( -#' formula_rhs = "caff_bin + age", -#' design = final_design, +#' formula_rhs = "sex + age", +#' design = design, #' var_to_test = "age", #' time_var = "stime", -#' status_var = "status" -#' ) -#' -#' # --- Example 2: Manual Safety Valve (if automatic extraction fails) --- -#' svycoxph_CE( -#' formula_rhs = "caff_bin + age", -#' design = final_design, -#' var_to_test = "caff_binYes", # Note: Check exact coef name -#' design_ids = ~psu, -#' design_weights = ~survey_weight, -#' design_strata = ~strata +#' status_var = "status", +#' n_intervals = 3, +#' print_main_model = FALSE, +#' print_split_summary = FALSE #' ) #' } diff --git a/R/svycoxph_CE_mi.R b/R/svycoxph_CE_mi.R index db96b55..376538e 100644 --- a/R/svycoxph_CE_mi.R +++ b/R/svycoxph_CE_mi.R @@ -17,9 +17,9 @@ #' @param tgroup_var A string for the time-interval variable created by #' \code{survSplit} (default is "tgroup"). #' @param time_var A string for the original 'stop' time variable in the -#' split data (default is "followup" or "stime"). -#' @param status_var A string for the event status variable -#' (default is "death" or "status"). +#' split data (default is "stime"). +#' @param status_var A string for the event status variable, coded 0/1 +#' (default is "status"). #' @param main_model_fit A pooled \code{svycoxph} fit object (class \code{mipo}). This is the #' main, un-split model. It is optional, but highly #' recommended, as the function will print its tidy @@ -56,7 +56,11 @@ #' summarized in the \strong{"Tidy Summary of All Time-Interval Models"} #' table. #' -#' @return A \code{ggplot} object showing the coefficient trend over time. +#' @return A \code{ggplot} object. The y-axis is the log hazard ratio for +#' \code{var_to_test} within each follow-up interval, pooled across imputations +#' (Rubin's rules) and plotted against time; a roughly flat trend supports the +#' proportional-hazards (constant-effect) assumption. This is a visual +#' diagnostic, not a formal hypothesis test. #' #' @importFrom mice pool #' @importFrom dplyr n @@ -65,147 +69,54 @@ #' @export #' #' @examples -#' \dontrun{ +#' \donttest{ +#' # MI version: test the constant-effect (PH) assumption across follow-up +#' # intervals, pooling over imputations with Rubin's rules. +#' if (requireNamespace("mice", quietly = TRUE) && +#' requireNamespace("mitools", quietly = TRUE)) { #' -#' # --- 1. Load Libraries --- -#' library(svyTable1) -#' library(dplyr) -#' library(survival) -#' library(survey) -#' library(mitools) -#' library(mice) -#' library(ggplot2) -#' library(knitr) +#' library(survival) # for Surv() / survSplit() in the formula below +#' data(nhanes_mortality, package = "svyTable1") +#' dat <- nhanes_mortality[nhanes_mortality$stime > 0, ] +#' dat$caff_bin <- factor(ifelse(dat$caff == "No consumption", "No", "Yes")) #' -#' # --- 2. Load data and create base design --- -#' data(nhanes_mortality) -#' analytic_design <- svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = nhanes_mortality, -#' nest = TRUE -#' ) +#' # Induce missingness, then impute (small m for a fast example only; for a +#' # real analysis use m >= 20 and check convergence). +#' set.seed(123) +#' dat$age[runif(nrow(dat)) < 0.10] <- NA +#' imp <- mice::mice(dat[, c("caff_bin", "sex", "age", "stime", "status", +#' "psu", "strata", "survey_weight")], +#' m = 2, maxit = 2, printFlag = FALSE, seed = 123) +#' long <- mice::complete(imp, "long", include = FALSE) #' -#' # --- 3. Create analysis data.frame and INDUCE MISSINGNESS in COVARIATES --- -#' set.seed(123) -#' data_with_miss <- analytic_design$variables %>% -#' filter(stime > 0) %>% -#' mutate( -#' # Create the exposure variable (complete) -#' caff_bin = factor( -#' ifelse(caff == "No consumption", "No", "Yes"), -#' levels = c("No", "Yes") -#' ), -#' # Induce 10% missingness in 'age' (a confounder) -#' age = ifelse(runif(n()) < 0.10, NA, age), -#' -#' # Induce 15% MAR missingness in 'bmi_cat' (a confounder) -#' bmi_cat = factor(ifelse(age > 50 & runif(n()) < 0.15, NA, as.character(bmi.cat))) +#' # Split follow-up at the median event time into two intervals. +#' cuts <- median(dat$stime[dat$status == 1], na.rm = TRUE) +#' long_split <- survSplit(Surv(stime, status) ~ ., +#' data = long, cut = cuts, +#' episode = "tgroup", id = "split_id") +#' design_split <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, nest = TRUE, +#' data = mitools::imputationList(split(long_split, long_split$.imp)) #' ) #' -#' print("--- Missingness Induced in Covariates ---") -#' print(mice::md.pattern(data_with_miss[, c("age", "bmi_cat", "caff_bin", "stime", "status")], -#' plot=FALSE)) -#' -#' # --- 4. Add Nelson-Aalen hazard (auxiliary variable for MICE) --- -#' data_with_miss$nelson_aalen <- nelsonaalen( -#' data_with_miss, -#' time = stime, -#' status = status -#' ) -#' -#' # --- 5. Run MICE --- -#' print("--- Starting MICE ---") -#' M_IMPUTATIONS <- 5 -#' MAX_ITERATIONS <- 5 -#' -#' pred_matrix <- make.predictorMatrix(data_with_miss) -#' vars_to_keep_as_is <- c("id", "survey.weight", "psu", "strata", -#' "stime", "status", "nelson_aalen", "caff_bin", "sex") -#' pred_matrix[, vars_to_keep_as_is] <- 0 -#' pred_matrix[vars_to_keep_as_is, ] <- 0 -#' -#' imputed_data <- mice( -#' data_with_miss, -#' m = M_IMPUTATIONS, -#' maxit = MAX_ITERATIONS, -#' predictorMatrix = pred_matrix, -#' method = 'pmm', -#' seed = 123, -#' printFlag = FALSE -#' ) -#' print("--- MICE Complete ---") -#' -#' # --- 6. Create Stacked Long-Format Data --- -#' impdata_long <- mice::complete(imputed_data, "long", include = FALSE) -#' -#' # --- 7. Fit the MAIN (Constant Effect) Model --- -#' print("--- Fitting Main Pooled (Constant Effect) Model ---") -#' -#' allImputations_main <- imputationList(split(impdata_long, impdata_long$.imp)) -#' design_main <- svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = allImputations_main, -#' nest = TRUE -#' ) -#' -#' my_formula <- "caff_bin + sex + age + bmi_cat" -#' main_formula <- as.formula(paste0("Surv(stime, status) ~ ", my_formula)) -#' -#' main_fit_list <- with(design_main, svycoxph(main_formula)) -#' main_fit_pooled <- pool(main_fit_list) -#' -#' # --- 8. Create the SPLIT-TIME Design --- -#' print("--- Creating Split-Time Design ---") -#' event_times <- data_with_miss$stime[data_with_miss$status == 1] -#' cuts <- quantile(event_times, probs = c(0.2, 0.4, 0.6, 0.8), na.rm = TRUE) -#' -#' impdata_long_split <- survSplit(Surv(stime, status) ~ ., -#' data = impdata_long, -#' cut = cuts, -#' episode = "tgroup", -#' id = "split_id") -#' -#' allImputations_split <- imputationList(split(impdata_long_split, -#' impdata_long_split$.imp)) -#' -#' design_split <- svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = allImputations_split, -#' nest = TRUE -#' ) -#' -#' # --- 9. Run the PH Test Function and Print the Plot --- -#' print("--- Calling svycoxph_CE_mi function to test 'caff_binYes' ---") -#' -#' my_ph_plot <- svycoxph_CE_mi( -#' formula_rhs = my_formula, -#' design_split = design_split, -#' var_to_test = "caff_binYes", # Test the exposure -#' tgroup_var = "tgroup", -#' time_var = "stime", -#' status_var = "status", -#' main_model_fit = main_fit_pooled, # Pass the main model here -#' print_split_summary = TRUE, -#' show_null_effect = TRUE -#' ) -#' -#' # Explicitly print the plot -#' print(my_ph_plot) -#' -#' } # End \dontrun{} +#' svycoxph_CE_mi( +#' formula_rhs = "caff_bin + sex + age", +#' design_split = design_split, +#' var_to_test = "caff_binYes", +#' tgroup_var = "tgroup", +#' time_var = "stime", +#' status_var = "status", +#' print_split_summary = FALSE +#' ) +#' } +#' } svycoxph_CE_mi <- function(formula_rhs, design_split, var_to_test, tgroup_var = "tgroup", - time_var = "followup", - status_var = "death", + time_var = "stime", + status_var = "status", main_model_fit = NULL, print_split_summary = TRUE, ...) { diff --git a/R/svydiag.R b/R/svydiag.R index 8ca7cdf..0ff0ae1 100644 --- a/R/svydiag.R +++ b/R/svydiag.R @@ -130,16 +130,16 @@ svydiag <- function(fit, p_threshold = 0.05, rse_threshold = 30) { ) %>% # Reorder and select the final columns for a clean output dplyr::select( - .data$Term, - .data$Estimate, - .data$SE, - .data$p.value, - .data$is_significant, - .data$CI_Lower, - .data$CI_Upper, - .data$CI_Width, - .data$RSE_percent, - .data$is_rse_high + "Term", + "Estimate", + "SE", + "p.value", + "is_significant", + "CI_Lower", + "CI_Upper", + "CI_Width", + "RSE_percent", + "is_rse_high" ) return(reliability_df) diff --git a/R/svygof.R b/R/svygof.R index 963cb86..aa24174 100644 --- a/R/svygof.R +++ b/R/svygof.R @@ -19,8 +19,10 @@ #' fitted probabilities. Defaults to 10 (deciles). #' #' @return -#' A `data.frame` containing the F-statistic, the numerator (df1) and -#' denominator (df2) degrees of freedom, and the p-value for the test. +#' A one-row `data.frame` with columns `F_statistic`, `df1` (numerator df), +#' `df2` (denominator df) and `p_value` for the Archer-Lemeshow design-based +#' goodness-of-fit test. A non-significant p-value (e.g. p > 0.05) is consistent +#' with adequate fit; it does not prove the model is correctly specified. #' #' @source #' The implementation is a formalized function based on the script and discussion @@ -32,86 +34,73 @@ #' @export #' #' @examples -#' \dontrun{ -#' # Ensure required packages are loaded -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE)) { +#' \donttest{ +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' # 1. Prepare Data -#' data(NHANESraw, package = "NHANES") -#' nhanes_data <- NHANESraw %>% -#' dplyr::filter(Age >= 20) %>% -#' dplyr::mutate(ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese"))) %>% -#' dplyr::filter(complete.cases(ObeseStatus, Age, Gender, Race1, -#' WTMEC2YR, SDMVPSU, SDMVSTRA)) +#' # The Archer-Lemeshow test uses a standard (non-replicate) survey design. +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # 2. Create a replicate design object -#' std_design <- survey::svydesign( -#' ids = ~SDMVPSU, -#' strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, -#' nest = TRUE, -#' data = nhanes_data -#' ) -#' rep_design <- survey::as.svrepdesign(std_design) +#' fit <- survey::svyglm(htn01 ~ age + sex + smoking, +#' design = design, family = quasibinomial()) #' -#' # 3. Fit a survey logistic regression model using the replicate design -#' fit_obesity_rep <- survey::svyglm( -#' ObeseStatus ~ Age + Gender + Race1, -#' design = rep_design, -#' family = quasibinomial() -#' ) -#' -#' # 4. Calculate the design-correct AUC -#' auc_results <- svyAUC(fit_obesity_rep, rep_design) -#' print(auc_results) -#' } +#' # A non-significant p-value is consistent with adequate fit (it does not +#' # prove the model is correctly specified). +#' svygof(fit, design, G = 10) #' } svygof <- function(fit, design, G = 10) { - # Get residuals and fitted values from the model + if (!inherits(fit, "svyglm")) { + stop("`fit` must be a survey logistic model of class 'svyglm'.", call. = FALSE) + } + if (!inherits(design, c("survey.design", "svyrep.design"))) { + stop("`design` must be a survey design ('survey.design' or 'svyrep.design').", + call. = FALSE) + } + + # Get response residuals and fitted probabilities from the model. resids <- stats::residuals(fit, type = "response") fitted_vals <- stats::fitted(fit) - # Create a data frame of model results, using row names to link back - model_data <- data.frame( - .id = names(resids), - r = resids, - f = fitted_vals - ) + # Align residuals/fitted to ALL rows of the design by row name, so that the + # decile grouping can be attached to the *existing* design object. Rows that + # were dropped while fitting (e.g. missing covariates) are left as NA and are + # excluded by na.omit when the grouped model is fitted. + rn <- rownames(design$variables) + n <- length(rn) + r_full <- rep(NA_real_, n) + f_full <- rep(NA_real_, n) + r_full[match(names(resids), rn)] <- as.numeric(resids) + f_full[match(names(fitted_vals), rn)] <- as.numeric(fitted_vals) - # Use the data directly from the design object, which is the most reliable source - data_with_res <- design$variables - data_with_res$.id <- rownames(data_with_res) - data_with_res <- merge(data_with_res, model_data, by = ".id", all.x = TRUE) - - # Create G groups based on fitted values - breaks <- stats::quantile(data_with_res$f, probs = seq(0, 1, 1 / G), na.rm = TRUE) + # Create G groups based on quantiles of the fitted probabilities. + breaks <- stats::quantile(f_full, probs = seq(0, 1, 1 / G), na.rm = TRUE) unique_breaks <- unique(breaks) - data_with_res$g <- cut(data_with_res$f, breaks = unique_breaks, include.lowest = TRUE) + if (length(unique_breaks) < 3) { + stop("Too few distinct fitted-probability groups for the goodness-of-fit ", + "test (the model's predictions are nearly constant). Try a smaller G.", + call. = FALSE) + } + g_full <- cut(f_full, breaks = unique_breaks, include.lowest = TRUE) - # Rebuild the design object using its internal components - new_design <- survey::svydesign( - ids = design$cluster, - strata = design$strata, - weights = design$weights, - data = data_with_res, - nest = isTRUE(design$nest) - ) + # Attach the residual and grouping variables to the EXISTING design via + # survey::update(). Unlike rebuilding the design from $cluster/$strata/$weights, + # this preserves finite population corrections, calibration, and replicate + # weights, so the Archer-Lemeshow test remains design-correct -- and it works + # for both 'survey.design' and 'svyrep.design' objects. + design <- stats::update(design, .gof_resid = r_full, .gof_group = g_full) - # Run the test - decile_model <- survey::svyglm(r ~ g, design = new_design, na.action = na.omit) - test_result <- survey::regTermTest(decile_model, ~g) + decile_model <- survey::svyglm(.gof_resid ~ .gof_group, design = design, + na.action = stats::na.omit) + test_result <- survey::regTermTest(decile_model, ~.gof_group) - # Return a tidy data frame - output <- data.frame( - F_statistic = test_result$Ftest[1], + data.frame( + F_statistic = as.numeric(test_result$Ftest[1]), df1 = test_result$df, df2 = test_result$ddf, - p_value = test_result$p + p_value = as.numeric(test_result$p) ) - - return(output) } diff --git a/R/svykmplot.R b/R/svykmplot.R index 6ac9d46..25b4821 100644 --- a/R/svykmplot.R +++ b/R/svykmplot.R @@ -17,6 +17,10 @@ #' use for the strata. #' @param show_pval Logical. If `TRUE`, calculates and displays the #' survey-weighted log-rank test p-value. +#' @param se Logical. If `TRUE` (default), computes survey-weighted confidence +#' bands via `survey::svykm(se = TRUE)` and draws them as shaded ribbons. The +#' standard-error calculation can be slow on large designs; set `se = FALSE` +#' to draw the survival curves only (much faster, no confidence bands). #' @param show_censor_marks Logical. If `TRUE`, displays censor markings (`+`) #' on the survival curves. #' @param base_font_size Base font size for the plot theme. @@ -48,54 +52,31 @@ #' @importFrom tools toTitleCase #' #' @examples -#' \dontrun{ -#' # This example uses the 'nhanes_mortality' dataset, +#' \donttest{ +#' library(survival) # for Surv() in the formula below +#' data(nhanes_mortality, package = "svyTable1") #' -#' if (requireNamespace("survey", quietly = TRUE)) { +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # 1. Load the data -#' data(nhanes_mortality) +#' # Survey-weighted Kaplan-Meier curves by caffeine consumption (women only). +#' design_female <- subset(design, sex == "Female") #' -#' # 2. Create the main survey design object -#' analytic_design <- survey::svydesign( -#' strata = ~strata, -#' id = ~psu, -#' weights = ~survey_weight, -#' data = nhanes_mortality, -#' nest = TRUE -#' ) +#' # se = FALSE draws curves only (fast); use se = TRUE to add confidence bands. +#' km <- svykmplot( +#' formula = Surv(stime, status) ~ caff, +#' design = design_female, +#' legend_title = "Caffeine consumption", +#' time_unit = "days", +#' time_breaks = seq(0, 240, by = 60), +#' show_pval = TRUE, +#' se = FALSE +#' ) #' -#' # 3. Create a subsetted design for females -#' design_female <- subset(analytic_design, sex == "Female") -#' -#' # 4. Define the formula -#' km_formula <- Surv(stime, status) ~ caff -#' -#' # 5. Define a 4-color palette -#' distinct_palette <- c("#377EB8", "#FF7F00", "#4DAF4A", "#E41A1C") -#' -#' # 6. Run the function -#' km_results_female <- svykmplot( -#' formula = km_formula, -#' design = design_female, -#' legend_title = "Caffeine Consumption", -#' time_unit = "days", # Use "days" so divisor is 1 -#' time_breaks = seq(0, 240, by = 60), -#' palette = distinct_palette, -#' show_pval = TRUE, -#' show_censor_marks = TRUE -#' ) -#' -#' # 7. Display the plot (and correct the x-axis label) -#' if (requireNamespace("ggplot2", quietly = TRUE)) { -#' km_results_female$plot + -#' ggplot2::labs(x = "Follow-up Time (Months)") -#' } -#' -#' # 8. Print the at-risk table -#' print(km_results_female$table) -#' -#' } +#' km$plot # a ggplot of the weighted survival curves +#' print(km$table) # the number-at-risk table #' } svykmplot <- function( formula, # e.g., Surv(stime, status) ~ age_meno @@ -106,6 +87,7 @@ svykmplot <- function( risk_table_title = "Number at Risk", palette = NULL, show_pval = TRUE, + se = TRUE, show_censor_marks = TRUE, base_font_size = 11, table_font_size = 3.5 @@ -157,7 +139,7 @@ svykmplot <- function( # --- 2. Run Survey Survival Calculations --- - svy_fit_object <- survey::svykm(formula, design = design_clean, se = TRUE) + svy_fit_object <- survey::svykm(formula, design = design_clean, se = se) n_strata <- length(svy_fit_object) strata_names <- names(svy_fit_object) @@ -177,7 +159,8 @@ svykmplot <- function( data.frame( time = svy_fit_object[[.x]]$time, surv = svy_fit_object[[.x]]$surv, - varlog = svy_fit_object[[.x]]$varlog, + # svykm only returns varlog when se = TRUE; use NA otherwise. + varlog = if (isTRUE(se)) svy_fit_object[[.x]]$varlog else NA_real_, strata = strata_names[.x] ) }) @@ -237,11 +220,18 @@ svykmplot <- function( curve_data, ggplot2::aes(x = .data$time, color = .data$strata, fill = .data$strata) ) + - ggplot2::geom_step(ggplot2::aes(y = .data$surv), lwd = 1) + - ggplot2::geom_ribbon( - ggplot2::aes(ymin = .data$lower_ci, ymax = .data$upper_ci), - alpha = 0.2, linetype = 0 - ) + + ggplot2::geom_step(ggplot2::aes(y = .data$surv), lwd = 1) + + # Confidence bands are only available (and only drawn) when se = TRUE. + if (isTRUE(se)) { + p_surv <- p_surv + + ggplot2::geom_ribbon( + ggplot2::aes(ymin = .data$lower_ci, ymax = .data$upper_ci), + alpha = 0.2, linetype = 0 + ) + } + + p_surv <- p_surv + ggplot2::theme_classic(base_size = base_font_size) + ggplot2::labs( y = "Survey-Weighted Survival Probability", @@ -334,7 +324,7 @@ svykmplot <- function( # --- 5. Combine and return --- printable_table <- table_data %>% - tidyr::pivot_wider(names_from = .data$time, values_from = .data$n.risk) %>% + tidyr::pivot_wider(names_from = "time", values_from = "n.risk") %>% dplyr::rename(!!legend_title := !!rlang::sym(strata_var)) return( diff --git a/R/svypooled.R b/R/svypooled.R index 721fb13..ccf144a 100644 --- a/R/svypooled.R +++ b/R/svypooled.R @@ -34,8 +34,12 @@ #' returns a table showing only the main exposure and lists covariates in a #' footnote. If `FALSE`, it returns a full table with all model terms. #' -#' @return An HTML table object of class `kableExtra`. This object can be -#' printed directly in R Markdown documents or saved. +#' @return A `kable` object (class `c("kableExtra", "knitr_kable")`): a character +#' vector of HTML with attributes, ready to print in an R Markdown report. In +#' "fallacy-safe" mode it contains one row per level of `main_exposure` (with the +#' adjustment variables named in a footnote); otherwise it contains all model +#' terms grouped by variable. Each cell shows the effect measure and its 95\% +#' confidence interval pooled across imputations using Rubin's rules. #' #' @importFrom dplyr select mutate case_when #' @importFrom stringr str_extract str_remove @@ -46,63 +50,32 @@ #' @export #' #' @examples -#' \dontrun{ -#' # Load required packages for the example -#' library(mice) -#' library(dplyr) -#' library(survey) -#' library(NHANES) +#' \donttest{ +#' data(nhanes_mortality, package = "svyTable1") +#' dat <- nhanes_mortality +#' dat$htn01 <- as.numeric(dat$htn == "Yes") +#' dat <- dat[, c("htn01", "sex", "age", "smoking", "psu", "strata", "survey_weight")] #' -#' # --- 1. Data Preparation --- -#' # We'll prepare a clean analytic dataset from the raw NHANES data. -#' # Note: We convert the outcome 'Obese' to a numeric 0/1 variable for svyglm. -#' data(NHANESraw, package = "NHANES") -#' nhanes_analytic <- NHANESraw %>% -#' dplyr::filter(Age >= 20 & WTMEC2YR > 0) %>% -#' mutate( -#' Obese_numeric = ifelse(BMI >= 30, 1, 0), -#' AgeCat = cut(Age, breaks = c(19, 39, 59, 80), labels = c("20-39", "40-59", "60-80")), -#' Smoke100 = factor(Smoke100, levels = c("No", "Yes")) -#' ) %>% -#' select(Obese_numeric, AgeCat, Smoke100, Education, SDMVPSU, SDMVSTRA, WTMEC2YR) +#' # Induce missingness, then multiply impute (small m for a fast example only; +#' # for a real analysis use m >= 20 and check convergence). +#' set.seed(1) +#' dat$age[sample(nrow(dat), 200)] <- NA +#' imp <- mice::mice(dat, m = 2, maxit = 2, printFlag = FALSE, seed = 1) #' -#' # --- 2. Perform Imputation and Pooled Analysis --- -#' # Set survey option to handle lonely PSUs, a common issue with NHANES data. -#' options(survey.lonely.psu = "adjust") +#' # Fit a survey-weighted model in each imputed dataset, then pool (Rubin's rules). +#' fits <- with(imp, survey::svyglm( +#' htn01 ~ sex + age + smoking, +#' design = survey::svydesign(id = ~psu, strata = ~strata, +#' weights = ~survey_weight, nest = TRUE, +#' data = data.frame(mget(ls()))), +#' family = quasibinomial() +#' )) +#' pooled <- mice::pool(fits) #' -#' # Impute the analytic dataset -#' imputed_data <- mice(nhanes_analytic, m = 2, maxit = 2, seed = 123, printFlag = FALSE) -#' -#' # Use with() to run a survey-weighted GLM on each imputed dataset -#' fit_imputed <- with(imputed_data, -#' svyglm(Obese_numeric ~ Smoke100 + AgeCat + Education, -#' design = svydesign(id = ~SDMVPSU, strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, nest = TRUE, data = .data), -#' family = quasibinomial()) -#' ) -#' -#' # Pool the results from all models into a single 'mipo' object -#' pooled_results <- pool(fit_imputed) -#' -#' # --- 3. Generate Tables with svypooled --- -#' # Example A: Create a "fallacy-safe" table (default) -#' svypooled( -#' pooled_model = pooled_results, -#' main_exposure = "Smoke100", -#' adj_var_names = c("AgeCat", "Education"), -#' measure = "OR", -#' title = "Adjusted Odds of Obesity (Fallacy-Safe)" -#' ) -#' -#' # Example B: Create a full table with all variables -#' svypooled( -#' pooled_model = pooled_results, -#' main_exposure = "Smoke100", -#' adj_var_names = c("AgeCat", "Education"), -#' measure = "OR", -#' title = "Adjusted Odds of Obesity (Full Table for Appendix)", -#' fallacy_safe = FALSE -#' ) +#' # Fallacy-safe table: shows only the exposure; covariates go in a footnote. +#' svypooled(pooled, main_exposure = "sex", +#' adj_var_names = c("age", "smoking"), measure = "OR", +#' title = "Adjusted odds of hypertension") #' } svypooled <- function(pooled_model, main_exposure, @@ -126,7 +99,7 @@ svypooled <- function(pooled_model, pattern <- paste(all_vars, collapse = "|") processed_results <- summary_df %>% - dplyr::select(.data$term, .data$estimate, .data$conf.low, .data$conf.high, .data$p.value) %>% + dplyr::select("term", "estimate", "conf.low", "conf.high", "p.value") %>% dplyr::mutate( group = stringr::str_extract(.data$term, pattern), Characteristic = stringr::str_remove(.data$term, pattern), @@ -136,7 +109,7 @@ svypooled <- function(pooled_model, TRUE ~ sprintf("%.3f", .data$p.value) ) ) %>% - dplyr::select(.data$group, .data$Characteristic, .data$Estimate_CI, .data$p_value_formatted) + dplyr::select("group", "Characteristic", "Estimate_CI", "p_value_formatted") # Conditionally filter for fallacy-safe output if (fallacy_safe) { diff --git a/R/svytable1.R b/R/svytable1.R index 9bac5af..3576d1f 100644 --- a/R/svytable1.R +++ b/R/svytable1.R @@ -35,8 +35,18 @@ #' reliability metrics. This includes the RSE and specific TRUE/FALSE columns #' (e.g., `fail_n_30`, `fail_eff_n_30`) for each suppression rule. Defaults to FALSE. #' -#' @return If `return_metrics` is FALSE (default), returns a data.frame. If TRUE, -#' returns a list with two elements: `formatted_table` and `reliability_metrics`. +#' @return If `return_metrics` is FALSE (default), a `data.frame` (the formatted +#' Table 1). Its first two columns are `Variable` and `Level`; the remaining +#' columns are `Overall` plus one column per level of `strata_var` (including a +#' `Missing` column when the stratifier has NAs). Cells hold unweighted counts +#' with weighted percentages for categorical variables and "mean (SD)" for +#' numeric variables; suppressed cells (when `reliability_checks = TRUE`) are +#' shown as `"*"`. If `return_metrics` is TRUE, a list with two elements: +#' `formatted_table` (the data.frame above) and `reliability_metrics` (a +#' data.frame with one row per evaluated cell giving `n`, `df`, `deff`, +#' `effective_n`, the confidence-interval bounds, `rse`, the logical +#' suppression flags `fail_n_30`/`fail_eff_n_30`/`fail_df_8`/`fail_ciw_30`/ +#' `fail_rciw_130`/`fail_rse_30`, and overall `suppressed`). #' #' @importFrom survey svytable svymean svyby svyvar svyciprop degf SE #' @import stats @@ -112,7 +122,11 @@ svytable1 <- function(design, strata_var, table_vars, } # --- 2. Table Generation --- - if (nrow(df) == 0) return(data.frame(Error = "Input data has 0 rows")) + if (nrow(df) == 0) { + stop("svytable1(): the survey design has 0 rows. The analytic dataset is ", + "empty - check your filter/drop_na/subset steps before building the design.", + call. = FALSE) + } if (any(is.na(df[[strata_var]]))) { strata_as_char <- as.character(df[[strata_var]]) @@ -229,7 +243,11 @@ svytable1 <- function(design, strata_var, table_vars, ci_low <- metrics$ci_low[level_index]; ci_high <- metrics$ci_high[level_index] pct_val <- metrics$prop[level_index]; se <- metrics$se[level_index] - effective_n <- if(!is.na(deff)) n / max(1, deff) else 0 + # Effective sample size = n / design effect. Guard only against a + # missing or non-positive deff (pathological in tiny cells); do NOT + # floor deff at 1, which would wrongly over-suppress cells that are + # more efficient than a simple random sample (deff < 1). + effective_n <- if (is.na(deff) || deff <= 0) 0 else n / deff ciw <- ci_high - ci_low rciw <- if(!is.na(pct_val) && pct_val > 0) (ciw / pct_val) * 100 else Inf rse <- if(!is.na(pct_val) && pct_val > 0) (se / pct_val) * 100 else Inf diff --git a/R/utils-pipe.R b/R/utils-pipe.R index 42de3c2..1d79f7f 100644 --- a/R/utils-pipe.R +++ b/R/utils-pipe.R @@ -7,6 +7,13 @@ #' @keywords internal #' @export #' @importFrom magrittr %>% +#' @return The result of calling the right-hand side with the left-hand side as +#' its first argument: \code{lhs \%>\% rhs} is equivalent to \code{rhs(lhs)}. +#' Used to chain a sequence of operations; the returned value and its class are +#' determined by the right-hand side expression. +#' @usage lhs \%>\% rhs +#' @param lhs A value to be piped into the right-hand side. +#' @param rhs A function call using the magrittr semantics. NULL utils::globalVariables(c( diff --git a/README.md b/README.md index 0ff5bc4..7018657 100644 --- a/README.md +++ b/README.md @@ -1,33 +1,30 @@ -# svyTable1: Create Publication-Ready Survey-Weighted Summary Tables +# svyTable1: Survey-Weighted Tables and Diagnostics for Epidemiology -**svyTable1** is an R package for creating publication-ready tables and visualizations from complex survey and epidemiologic analyses. It streamlines the creation of descriptive “Table 1” summaries, supports multiply imputed regression models, provides diagnostic and goodness-of-fit tools, enables interaction effect reporting, and supports survey-weighted survival analysis. +**svyTable1** is a teaching/helper R package for the UBC SPPH 604 course and the +[EpiMethods](https://ehsanx.github.io/EpiMethods/) book. It helps students build +publication-ready tables and diagnostics from complex survey data (such as +NHANES): descriptive "Table 1" summaries, survey-weighted regression and survival +diagnostics, additive-interaction reporting, and pooled multiply-imputed models. + +All examples below use the bundled `nhanes_mortality` dataset, so they run without +any additional downloads. --- ## ✨ Key Features -- **Survey Design Integration:** Built on `svydesign` objects from the **survey** package. -- **Automatic Missing Data Handling:** Automatically reports missing values for each variable. -- **Best Practice Reporting:** Defaults to “mixed mode” with unweighted Ns and weighted estimates. -- **Reliability Checks:** Implements NCHS Data Presentation Standards for Proportions. -- **Flexible Output Modes:** `"mixed"`, `"weighted"`, and `"unweighted"`. -- **Survival Plotting:** `svykmplot()` for publication-ready survey-weighted Kaplan Meier plots with number at risk. -- **Regression Diagnostics:** `svydiag()` to evaluate coefficient reliability. -- **Goodness-of-Fit:** `svygof()` to perform Archer-Lemeshow GOF tests for survey logistic regression. -- **Design-Correct AUC:** `svyAUC()` for valid AUC estimation under complex survey designs. -- **Multiple Imputation Support:** `svypooled()` generates fallacy-safe pooled regression tables from `mice`. -- **Interaction Analysis:** - - *Joint effects*: `jointeffects` - - *Simple effects*: `inteffects` - - *Additive interaction*: `addint`, `addintlist` -- **Survey Cox Regression with Contrast Estimation:** - `svycoxph_CE()` produces hazard ratio contrasts for groups or combinations. -- **MI-Compatible Cox Regression with Contrasts:** - `svycoxph_CE_mi()` extends the above to multiply imputed datasets. -- **Plotting Interaction Effects:** - `plotint()` generates publication-ready plots for interaction effects estimated with `addint`, `jointeffects`, or `inteffects`. -- **Automated Interaction Reporting:** - `reportint()` creates narrative-ready summaries of interaction effects, RERI, AP, and joint effects. +- **Built on the `survey` package** — works with `svydesign` objects. +- **Descriptive tables** — `svytable1()` with unweighted *n* + weighted %, automatic + missing-data handling, and optional NCHS reliability suppression. +- **Regression diagnostics** — `svydiag()` (coefficient reliability), `svygof()` + (Archer–Lemeshow goodness-of-fit), `svyAUC()` (design-correct AUC). +- **Interaction analysis** — additive interaction (`addint()`, `addintlist()`), + joint effects (`jointeffects()`), simple effects (`inteffects()`), a multi-panel + report (`reportint()`), and interaction plots (`plotint()`). +- **Survival analysis** — survey-weighted Kaplan–Meier plots with number-at-risk + (`svykmplot()`) and constant-effect / proportional-hazards diagnostics + (`svycoxph_CE()`, `svycoxph_CE_mi()`). +- **Multiple imputation** — fallacy-safe pooled tables from `mice` (`svypooled()`). --- @@ -40,117 +37,150 @@ devtools::install_github("ehsanx/svyTable1", build_vignettes = TRUE, dependencie --- -## 🔍 Usage Examples (Selected) - -Below are a few key examples. See vignettes for complete coverage. - ---- - -## Example: Interaction Reporting with reportint() +## 🚀 Quick start ```r library(svyTable1) library(survey) -library(dplyr) -data(NHANESraw) +options(survey.lonely.psu = "adjust") # correct SEs for NHANES-style designs -analytic <- NHANESraw %>% - filter(Age >= 20) %>% - mutate( - Obese = factor(ifelse(BMI >= 30, "Yes", "No")), - Hypertension = factor(ifelse(BPSysAve >= 130, "Yes", "No")) - ) %>% - drop_na(Age, Race1, Obese, Hypertension, SDMVPSU, SDMVSTRA, WTMEC2YR) +data(nhanes_mortality) +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") design <- svydesign( - id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, - nest = TRUE, data = analytic + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality ) +``` + +### Descriptive "Table 1" -model_int <- svyglm( - Hypertension ~ Obese * Race1 + Age, - design = design, - family = quasibinomial() +```r +svytable1( + design = design, + strata_var = "htn", + table_vars = c("age", "sex", "smoking", "bmi.cat") ) -reportint(model_int, factor1_name = "Obese", factor2_name = "Race1") +# With NCHS reliability checks and a detailed metrics table: +svytable1( + design = design, + strata_var = "htn", + table_vars = c("age", "sex", "smoking", "bmi.cat"), + reliability_checks = TRUE, + return_metrics = TRUE +) ``` --- -## Example: Plotting Interaction Effects with plotint() +## 🔍 Regression diagnostics ```r -interaction_results <- addint( - model = model_int, - factor1_name = "Obese", - factor2_name = "Race1", - measures = "all" -) +fit <- svyglm(htn01 ~ age + sex + smoking, design = design, family = quasibinomial()) + +# Per-coefficient reliability (SE, p, CI width, RSE) +svydiag(fit) -plotint(interaction_results, measure = "RERI") +# Archer-Lemeshow goodness-of-fit (use a standard, non-replicate design) +svygof(fit, design) + +# Design-correct AUC (needs a replicate-weights design); CI is on the logit scale +rep_design <- as.svrepdesign(design) +fit_rep <- svyglm(htn01 ~ age + sex + smoking, design = rep_design, + family = quasibinomial()) +svyAUC(fit_rep, rep_design) ``` --- -## Example: Survey Cox Regression with Contrast Estimation +## 🔗 Interaction analysis ```r -library(survival) +int_fit <- svyglm(htn01 ~ sex * insulin + age, design = design, + family = quasibinomial()) -cox_model <- svycoxph( - Surv(followup_time, death) ~ Obese + Age + Race1, - design = design -) +# Additive interaction measures (RERI, AP, Synergy index) for every level pair +addintlist(int_fit, factor1_name = "sex", factor2_name = "insulin", measures = "all") -svycoxph_CE( - model = cox_model, - contrast_var = "Obese", - ref = "No", - target = "Yes" -) +# Joint effects of the two factors +jointeffects(int_fit, "sex", "insulin") + +# Multi-panel interaction report (joint, stratum-specific, additive, multiplicative) +report <- reportint(int_fit, factor1_name = "sex", factor2_name = "insulin", + output = "list") +report$Interaction_report ``` ---- +> **Note on scale.** `addint()`/`addintlist()` compute RERI/AP/S by exponentiating +> the model coefficients. For logistic models these are odds ratios, which only +> approximate risk ratios when the outcome is rare. Interpret additive-interaction +> measures on the risk-ratio scale with care for common outcomes. -## Example: Cox Contrasts with Multiple Imputation +Plotting an interaction with a continuous moderator: ```r -library(mice) - -imp <- mice(analytic, m = 5, seed = 123, printFlag = FALSE) +int_fit2 <- svyglm(htn01 ~ insulin * age, design = design, family = quasibinomial()) -fit_list <- with(imp, { - d <- svydesign(id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, - nest = TRUE, data = complete(data)) - svycoxph(Surv(followup_time, death) ~ Obese + Age + Race1, design = d) -}) - -svycoxph_CE_mi(fit_list, contrast_var = "Obese", ref = "No", target = "Yes") +plotint(model = int_fit2, effect = "insulin", moderator = "age", + data = nhanes_mortality) ``` --- -## Example: Creating Interaction Tables +## ⏳ Survival analysis ```r -addintlist( - model = model_int, - factor1_name = "Obese", - factor2_name = "Race1", - measures = "all" +library(survival) + +# Survey-weighted Kaplan-Meier curves with a number-at-risk table. +# se = TRUE adds confidence bands but is slow on large designs; use se = FALSE +# to draw curves only. +km <- svykmplot( + formula = Surv(stime, status) ~ caff, + design = subset(design, sex == "Female"), + time_unit = "days", + time_breaks = seq(0, 240, by = 60), + show_pval = TRUE, + se = FALSE +) +km$plot +print(km$table) + +# Constant-effect / proportional-hazards diagnostic for a Cox model +svycoxph_CE( + formula_rhs = "sex + age", + design = design, + var_to_test = "age", + time_var = "stime", + status_var = "status", + n_intervals = 3 ) ``` +A multiple-imputation version of the Cox diagnostic +(`svycoxph_CE_mi()`) and a fallacy-safe pooled table (`svypooled()`) are +demonstrated in `?svycoxph_CE_mi`, `?svypooled`, and the package vignettes. + --- -## 🤝 Contributing +## 📚 Learn more + +```r +browseVignettes("svyTable1") +``` -Contributions are welcome. +The vignettes follow the SPPH 604 / EpiMethods workflow: descriptive tables → +model diagnostics → interaction analysis → survival analysis → multiple imputation. --- +## 🤝 Contributing + +Contributions and bug reports are welcome via the +[issue tracker](https://github.com/ehsanx/svyTable1/issues). + ## 📜 License MIT License. diff --git a/man/addint.Rd b/man/addint.Rd index 821f1ff..7ef6725 100644 --- a/man/addint.Rd +++ b/man/addint.Rd @@ -61,97 +61,32 @@ Confidence intervals for S are calculated on the log scale and then exponentiate } \examples{ \donttest{ -# --- Load required libraries for the example --- -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE) && - requireNamespace("tidyr", quietly = TRUE) && - requireNamespace("msm", quietly = TRUE)) { +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - library(survey) - library(NHANES) - library(dplyr) - library(tidyr) - library(msm) - - # --- 1. Data Preparation (NHANES Example) --- - data(NHANESraw) - - vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", - "SDMVPSU", "SDMVSTRA", "WTMEC2YR") - - nhanes_adults_processed <- NHANESraw \%>\% - filter(Age >= 20) \%>\% - mutate( - ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese")), - Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), - levels = c("No", "Yes")), - Race1 = relevel(as.factor(Race1), ref = "White") - ) \%>\% - select(all_of(vars_needed)) \%>\% - drop_na() - - adult_design_binary <- svydesign( - id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, - nest = TRUE, data = nhanes_adults_processed - ) - - # --- 2a. Fit Joint Variable Model --- - adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE) - ) - adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese") - ) - - joint_model_logit <- svyglm( - Hypertension_130 ~ Race1_ObeseStatus + Age, - design = adult_design_binary, family = quasibinomial() - ) - - # --- 2b. Fit Interaction Term Model --- - interaction_model_logit <- svyglm( - Hypertension_130 ~ Race1 * ObeseStatus + Age, - design = adult_design_binary, family = quasibinomial() - ) - - # --- 3. Calculate Additive Interaction (Black vs White * Obese vs Not Obese) --- - - # --- Using the Joint Model --- - joint_coef_names_black_obese <- list( - exp1_level = "Race1_ObeseStatusBlack_Not Obese", # A=1, B=0 - exp2_level = "Race1_ObeseStatusWhite_Obese", # A=0, B=1 - both_levels = "Race1_ObeseStatusBlack_Obese" # A=1, B=1 - ) - - results_joint <- addint( - model = joint_model_logit, - type = "joint", - coef_names = joint_coef_names_black_obese, - measures = "all" - ) - print("--- Results from Joint Model (Black vs Obese Interaction) ---") - print(results_joint) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # --- Using the Interaction Model --- - interaction_coef_names_black_obese <- list( - exp1_coef = "Race1Black", # Main effect A - exp2_coef = "ObeseStatusObese", # Main effect B - inter_coef = "Race1Black:ObeseStatusObese" # Interaction A:B - ) +# Interaction-term model: additive interaction of sex and insulin use. +fit <- survey::svyglm(htn01 ~ sex * insulin + age, + design = design, family = quasibinomial()) - results_interact <- addint( - model = interaction_model_logit, - type = "interaction", - coef_names = interaction_coef_names_black_obese, - measures = "all" - ) - print("--- Results from Interaction Model (Black vs Obese Interaction) ---") - print(results_interact) +# Coefficient names must match names(coef(fit)) exactly. +addint( + model = fit, + type = "interaction", + coef_names = list( + exp1_coef = "sexFemale", + exp2_coef = "insulinYes", + inter_coef = "sexFemale:insulinYes" + ), + measures = "all" +) -} else { - print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -} +# RERI/AP/S can also be computed from a joint-variable model with +# type = "joint" (see the 'coef_names' argument). For logistic models these +# use odds ratios, which approximate risk ratios only when the outcome is rare. } } diff --git a/man/addintlist.Rd b/man/addintlist.Rd index 8a7fbcb..5e33ef1 100644 --- a/man/addintlist.Rd +++ b/man/addintlist.Rd @@ -10,7 +10,8 @@ addintlist( factor2_name, measures = "all", conf.level = 0.95, - digits = 3 + digits = 3, + verbose = FALSE ) } \arguments{ @@ -30,6 +31,9 @@ Options: `"RERI"`, `"AP"`, `"S"`, or `"all"`. Default is `"all"`.} \item{digits}{Integer. Number of decimal places for rounding estimates and CIs in the final table. Default is 3. (Set to NULL for no rounding).} + +\item{verbose}{Logical. If `TRUE`, prints a progress message for each +factor-level combination as it is calculated. Defaults to `FALSE`.} } \value{ A `tibble` (data frame) summarizing the additive interaction results. diff --git a/man/inteffects.Rd b/man/inteffects.Rd index 960e266..aad984e 100644 --- a/man/inteffects.Rd +++ b/man/inteffects.Rd @@ -69,86 +69,32 @@ with consistent column structure, including scale-appropriate standard errors. } \examples{ \donttest{ -# --- Load required libraries for the example --- -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE) && - requireNamespace("tidyr", quietly = TRUE) && - requireNamespace("msm", quietly = TRUE)) { # msm needed for SE calculation +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - library(survey) - library(NHANES) - library(dplyr) - library(tidyr) - library(msm) # Load msm - - # --- 1. Data Preparation (NHANES Example) --- - data(NHANESraw) - - vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", - "SDMVPSU", "SDMVSTRA", "WTMEC2YR") - - nhanes_adults_processed <- NHANESraw \%>\% - filter(Age >= 20) \%>\% - mutate( - ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese")), - Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), - levels = c("No", "Yes")), - Race1 = relevel(as.factor(Race1), ref = "White") - ) \%>\% - select(all_of(vars_needed)) \%>\% - drop_na() - - adult_design_binary <- svydesign( - id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, - nest = TRUE, data = nhanes_adults_processed - ) - - # --- 2. Create Joint Variable and Fit Joint Model --- - adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE) - ) - adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese") - ) - - joint_model_logit <- svyglm( - Hypertension_130 ~ Race1_ObeseStatus + Age, - design = adult_design_binary, family = quasibinomial() - ) +# Build a joint variable from two factors, then fit a joint-variable model. +nhanes_mortality$sex_insulin <- interaction( + nhanes_mortality$sex, nhanes_mortality$insulin, sep = "_" +) - # --- 3. Run the function --- - f1_levels <- levels(adult_design_binary$variables$Race1) - f2_levels <- levels(adult_design_binary$variables$ObeseStatus) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # --- Example 1: Output on Ratio scale --- - simple_effects_ratio <- inteffects( - joint_model = joint_model_logit, - joint_var_name = "Race1_ObeseStatus", - factor1_name = "Race1", - factor2_name = "ObeseStatus", - factor1_levels = f1_levels, - factor2_levels = f2_levels, - level_separator = "_", - scale = "ratio" - ) - print("--- Output on Ratio Scale ---") - print(simple_effects_ratio, n = 50) +joint_fit <- survey::svyglm(htn01 ~ sex_insulin + age, + design = design, family = quasibinomial()) - # --- Example 2: Output on Log scale --- - simple_effects_log <- inteffects( - joint_model = joint_model_logit, - joint_var_name = "Race1_ObeseStatus", - factor1_name = "Race1", - factor2_name = "ObeseStatus", - factor1_levels = f1_levels, - factor2_levels = f2_levels, - level_separator = "_", - scale = "log" - ) - print("--- Output on Log Scale ---") - print(simple_effects_log, n = 50) -} +# Stratum-specific (simple) effects reconstructed from the joint model. +inteffects( + joint_model = joint_fit, + joint_var_name = "sex_insulin", + factor1_name = "sex", + factor2_name = "insulin", + factor1_levels = c("Male", "Female"), # reference level first + factor2_levels = c("No", "Yes"), # reference level first + level_separator = "_", + scale = "ratio" +) } } diff --git a/man/pipe.Rd b/man/pipe.Rd index c5480cc..5ec80d6 100644 --- a/man/pipe.Rd +++ b/man/pipe.Rd @@ -3,6 +3,20 @@ \name{\%>\%} \alias{\%>\%} \title{Pipe operator} +\usage{ +lhs \%>\% rhs +} +\arguments{ +\item{lhs}{A value to be piped into the right-hand side.} + +\item{rhs}{A function call using the magrittr semantics.} +} +\value{ +The result of calling the right-hand side with the left-hand side as + its first argument: \code{lhs \%>\% rhs} is equivalent to \code{rhs(lhs)}. + Used to chain a sequence of operations; the returned value and its class are + determined by the right-hand side expression. +} \description{ See \code{magrittr::\link[magrittr]{\%>\%}} for details. } diff --git a/man/reportint.Rd b/man/reportint.Rd index 5c9a5c1..3c609b9 100644 --- a/man/reportint.Rd +++ b/man/reportint.Rd @@ -69,70 +69,27 @@ It takes a **single pre-fitted model** as input and works for `glm`, This function requires the 'svyTable1' and 'Publish' packages. } \examples{ -\dontrun{ -# --- Load required libraries for the example --- -if (requireNamespace("svyTable1", quietly = TRUE) && - requireNamespace("Publish", quietly = TRUE) && - requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE) && - requireNamespace("tidyr", quietly = TRUE) && - requireNamespace("broom", quietly = TRUE)) { +\donttest{ +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - library(svyTable1) - library(Publish) - library(survey) - library(NHANES) - library(dplyr) - library(tidyr) - library(broom) - - # --- 1. Data Preparation (NHANES Example) --- - data(NHANESraw) - - vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130", - "SDMVPSU", "SDMVSTRA", "WTMEC2YR") - - nhanes_adults_processed <- NHANESraw \%>\% - dplyr::filter(Age >= 20) \%>\% - dplyr::mutate( - ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese")), - Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"), - levels = c("No", "Yes")), - Race1 = relevel(as.factor(Race1), ref = "White") - ) \%>\% - dplyr::select(dplyr::all_of(vars_needed)) \%>\% - tidyr::drop_na() - - adult_design_binary <- svydesign( - id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, - nest = TRUE, data = nhanes_adults_processed - ) - - # --- 2. FIT ONLY THE INTERACTION MODEL --- - - interaction_model_fit <- svyglm( - Hypertension_130 ~ Race1 * ObeseStatus + Age, - design = adult_design_binary, family = quasibinomial() - ) - - # --- 3. Run the wrapper function --- - - report_list <- reportint( - interaction_model = interaction_model_fit - ) - - # You can then print the individual panel tables: - # print(report_list$joint_effects) - # print(report_list$stratum_specific_effects) - - # --- OR --- +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # Print the single, publication-ready summary tables: - # print(report_list$effect_modification_report) - # print(report_list$Interaction_report) +# Fit a model with a two-factor interaction (sex x insulin). +interaction_fit <- survey::svyglm( + htn01 ~ sex * insulin + age, + design = design, family = quasibinomial() +) -} +# The default output = "list" returns the panel data frames (no rendering). +report <- reportint(interaction_fit, + factor1_name = "sex", factor2_name = "insulin", + output = "list") +report$joint_effects +report$additive_interaction +report$Interaction_report } } diff --git a/man/svyAUC.Rd b/man/svyAUC.Rd index a5cbc3e..c22904f 100644 --- a/man/svyAUC.Rd +++ b/man/svyAUC.Rd @@ -14,9 +14,13 @@ svyAUC(fit, design, plot = FALSE) \item{plot}{A logical value. If `TRUE`, an ROC curve is plotted. Defaults to `FALSE`.} } \value{ -If `plot = FALSE` (default), returns a `data.frame` with the AUC, SE, and 95% CI. -If `plot = TRUE`, invisibly returns a list containing the summary `data.frame` -and another `data.frame` with the ROC curve coordinates (TPR and FPR). +If `plot = FALSE` (default), a one-row `data.frame` with columns `AUC`, `SE`, +`CI_Lower` and `CI_Upper`. The standard error is computed across the design's +replicate weights and the 95\% confidence interval is built on the logit scale +(so it stays within `[0, 1]`); both are therefore conditional on the +replication scheme used to build `design`. If `plot = TRUE`, an ROC curve is +drawn and a list is returned invisibly with the summary `data.frame` and a +`data.frame` of ROC coordinates (TPR and FPR). } \description{ This function calculates the Area Under the Curve (AUC) and its design-correct @@ -24,41 +28,21 @@ standard error and 95% confidence interval. It can also generate a plot of the weighted ROC curve. } \examples{ -\dontrun{ -# Ensure required packages are loaded -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE)) { +\donttest{ +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - # 1. Prepare Data - data(NHANESraw, package = "NHANES") - nhanes_data <- NHANESraw \%>\% - dplyr::filter(Age >= 20) \%>\% - dplyr::mutate(ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese"))) \%>\% - dplyr::filter(complete.cases(ObeseStatus, Age, Gender, Race1, - WTMEC2YR, SDMVPSU, SDMVSTRA)) +# svyAUC requires a replicate-weights design. +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) +rep_design <- survey::as.svrepdesign(design) - # 2. Create a replicate design object - std_design <- survey::svydesign( - ids = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_data - ) - rep_design <- survey::as.svrepdesign(std_design) +fit <- survey::svyglm(htn01 ~ age + sex + smoking, + design = rep_design, family = quasibinomial()) - # 3. Fit a survey logistic regression model using the replicate design - fit_obesity_rep <- survey::svyglm( - ObeseStatus ~ Age + Gender + Race1, - design = rep_design, - family = quasibinomial() - ) - - # 4. Calculate the design-correct AUC - auc_results <- svyAUC(fit_obesity_rep, rep_design) - print(auc_results) -} +# AUC with a design-correct SE and a logit-scale 95\% CI. +svyAUC(fit, rep_design) } } diff --git a/man/svycoxph_CE.Rd b/man/svycoxph_CE.Rd index a3098c2..69de07d 100644 --- a/man/svycoxph_CE.Rd +++ b/man/svycoxph_CE.Rd @@ -71,7 +71,10 @@ coefficients for every time interval.} \item{show_null_effect}{Logical. If TRUE, draws a dashed red line at y=0.} } \value{ -A \code{ggplot} object showing the coefficient trend over time. +A \code{ggplot} object. The y-axis is the estimated log hazard ratio + for \code{var_to_test} within each follow-up interval, plotted against time; + a roughly flat trend supports the proportional-hazards (constant-effect) + assumption. This is a visual diagnostic, not a formal hypothesis test. } \description{ This function provides a valid alternative to \code{cox.zph()} for complex survey @@ -95,52 +98,25 @@ before splitting. If this fails, the user can manually provide the \code{design_ids}, \code{design_weights}, and \code{design_strata} arguments. } \examples{ -\dontrun{ -# --- 1. Load data and create a design --- -library(survey) -library(dplyr) +\donttest{ data(nhanes_mortality, package = "svyTable1") -# Create a base design -analytic_design <- svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = nhanes_mortality, - nest = TRUE +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality ) -# Prepare data (filter to >0 time) -data_clean <- analytic_design$variables \%>\% - dplyr::filter(stime > 0) \%>\% - mutate(caff_bin = ifelse(caff == "No consumption", "No", "Yes")) - -# Create final design object -final_design <- svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = data_clean, - nest = TRUE -) - -# --- Example 1: Standard Usage --- +# Visual check of the constant-effect (proportional hazards) assumption for +# 'age': a roughly flat trend across follow-up intervals supports it. svycoxph_CE( - formula_rhs = "caff_bin + age", - design = final_design, + formula_rhs = "sex + age", + design = design, var_to_test = "age", time_var = "stime", - status_var = "status" -) - -# --- Example 2: Manual Safety Valve (if automatic extraction fails) --- -svycoxph_CE( - formula_rhs = "caff_bin + age", - design = final_design, - var_to_test = "caff_binYes", # Note: Check exact coef name - design_ids = ~psu, - design_weights = ~survey_weight, - design_strata = ~strata + status_var = "status", + n_intervals = 3, + print_main_model = FALSE, + print_split_summary = FALSE ) } } diff --git a/man/svycoxph_CE_mi.Rd b/man/svycoxph_CE_mi.Rd index ed288b5..6c098c9 100644 --- a/man/svycoxph_CE_mi.Rd +++ b/man/svycoxph_CE_mi.Rd @@ -9,8 +9,8 @@ svycoxph_CE_mi( design_split, var_to_test, tgroup_var = "tgroup", - time_var = "followup", - status_var = "death", + time_var = "stime", + status_var = "status", main_model_fit = NULL, print_split_summary = TRUE, ... @@ -29,10 +29,10 @@ has *already been split* using \code{survSplit()}.} \code{survSplit} (default is "tgroup").} \item{time_var}{A string for the original 'stop' time variable in the -split data (default is "followup" or "stime").} +split data (default is "stime").} -\item{status_var}{A string for the event status variable -(default is "death" or "status").} +\item{status_var}{A string for the event status variable, coded 0/1 +(default is "status").} \item{main_model_fit}{A pooled \code{svycoxph} fit object (class \code{mipo}). This is the main, un-split model. It is optional, but highly @@ -48,7 +48,11 @@ of all pooled coefficients from all interval models.} \code{show_null_effect} (TRUE/FALSE).} } \value{ -A \code{ggplot} object showing the coefficient trend over time. +A \code{ggplot} object. The y-axis is the log hazard ratio for + \code{var_to_test} within each follow-up interval, pooled across imputations + (Rubin's rules) and plotted against time; a roughly flat trend supports the + proportional-hazards (constant-effect) assumption. This is a visual + diagnostic, not a formal hypothesis test. } \description{ This function provides a valid alternative to \code{cox.zph()} for complex survey @@ -84,138 +88,45 @@ This function helps you decide which model to report in your paper. table. } \examples{ -\dontrun{ - -# --- 1. Load Libraries --- -library(svyTable1) -library(dplyr) -library(survival) -library(survey) -library(mitools) -library(mice) -library(ggplot2) -library(knitr) - -# --- 2. Load data and create base design --- -data(nhanes_mortality) -analytic_design <- svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = nhanes_mortality, - nest = TRUE -) - -# --- 3. Create analysis data.frame and INDUCE MISSINGNESS in COVARIATES --- -set.seed(123) -data_with_miss <- analytic_design$variables \%>\% - filter(stime > 0) \%>\% - mutate( - # Create the exposure variable (complete) - caff_bin = factor( - ifelse(caff == "No consumption", "No", "Yes"), - levels = c("No", "Yes") - ), - # Induce 10\% missingness in 'age' (a confounder) - age = ifelse(runif(n()) < 0.10, NA, age), - - # Induce 15\% MAR missingness in 'bmi_cat' (a confounder) - bmi_cat = factor(ifelse(age > 50 & runif(n()) < 0.15, NA, as.character(bmi.cat))) +\donttest{ +# MI version: test the constant-effect (PH) assumption across follow-up +# intervals, pooling over imputations with Rubin's rules. +if (requireNamespace("mice", quietly = TRUE) && + requireNamespace("mitools", quietly = TRUE)) { + + library(survival) # for Surv() / survSplit() in the formula below + data(nhanes_mortality, package = "svyTable1") + dat <- nhanes_mortality[nhanes_mortality$stime > 0, ] + dat$caff_bin <- factor(ifelse(dat$caff == "No consumption", "No", "Yes")) + + # Induce missingness, then impute (small m for a fast example only; for a + # real analysis use m >= 20 and check convergence). + set.seed(123) + dat$age[runif(nrow(dat)) < 0.10] <- NA + imp <- mice::mice(dat[, c("caff_bin", "sex", "age", "stime", "status", + "psu", "strata", "survey_weight")], + m = 2, maxit = 2, printFlag = FALSE, seed = 123) + long <- mice::complete(imp, "long", include = FALSE) + + # Split follow-up at the median event time into two intervals. + cuts <- median(dat$stime[dat$status == 1], na.rm = TRUE) + long_split <- survSplit(Surv(stime, status) ~ ., + data = long, cut = cuts, + episode = "tgroup", id = "split_id") + design_split <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, nest = TRUE, + data = mitools::imputationList(split(long_split, long_split$.imp)) ) -print("--- Missingness Induced in Covariates ---") -print(mice::md.pattern(data_with_miss[, c("age", "bmi_cat", "caff_bin", "stime", "status")], - plot=FALSE)) - -# --- 4. Add Nelson-Aalen hazard (auxiliary variable for MICE) --- -data_with_miss$nelson_aalen <- nelsonaalen( - data_with_miss, - time = stime, - status = status -) - -# --- 5. Run MICE --- -print("--- Starting MICE ---") -M_IMPUTATIONS <- 5 -MAX_ITERATIONS <- 5 - -pred_matrix <- make.predictorMatrix(data_with_miss) -vars_to_keep_as_is <- c("id", "survey.weight", "psu", "strata", - "stime", "status", "nelson_aalen", "caff_bin", "sex") -pred_matrix[, vars_to_keep_as_is] <- 0 -pred_matrix[vars_to_keep_as_is, ] <- 0 - -imputed_data <- mice( - data_with_miss, - m = M_IMPUTATIONS, - maxit = MAX_ITERATIONS, - predictorMatrix = pred_matrix, - method = 'pmm', - seed = 123, - printFlag = FALSE -) -print("--- MICE Complete ---") - -# --- 6. Create Stacked Long-Format Data --- -impdata_long <- mice::complete(imputed_data, "long", include = FALSE) - -# --- 7. Fit the MAIN (Constant Effect) Model --- -print("--- Fitting Main Pooled (Constant Effect) Model ---") - -allImputations_main <- imputationList(split(impdata_long, impdata_long$.imp)) -design_main <- svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = allImputations_main, - nest = TRUE -) - -my_formula <- "caff_bin + sex + age + bmi_cat" -main_formula <- as.formula(paste0("Surv(stime, status) ~ ", my_formula)) - -main_fit_list <- with(design_main, svycoxph(main_formula)) -main_fit_pooled <- pool(main_fit_list) - -# --- 8. Create the SPLIT-TIME Design --- -print("--- Creating Split-Time Design ---") -event_times <- data_with_miss$stime[data_with_miss$status == 1] -cuts <- quantile(event_times, probs = c(0.2, 0.4, 0.6, 0.8), na.rm = TRUE) - -impdata_long_split <- survSplit(Surv(stime, status) ~ ., - data = impdata_long, - cut = cuts, - episode = "tgroup", - id = "split_id") - -allImputations_split <- imputationList(split(impdata_long_split, - impdata_long_split$.imp)) - -design_split <- svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = allImputations_split, - nest = TRUE -) - -# --- 9. Run the PH Test Function and Print the Plot --- -print("--- Calling svycoxph_CE_mi function to test 'caff_binYes' ---") - -my_ph_plot <- svycoxph_CE_mi( - formula_rhs = my_formula, - design_split = design_split, - var_to_test = "caff_binYes", # Test the exposure - tgroup_var = "tgroup", - time_var = "stime", - status_var = "status", - main_model_fit = main_fit_pooled, # Pass the main model here - print_split_summary = TRUE, - show_null_effect = TRUE -) - -# Explicitly print the plot -print(my_ph_plot) - -} # End \dontrun{} + svycoxph_CE_mi( + formula_rhs = "caff_bin + sex + age", + design_split = design_split, + var_to_test = "caff_binYes", + tgroup_var = "tgroup", + time_var = "stime", + status_var = "status", + print_split_summary = FALSE + ) +} +} } diff --git a/man/svygof.Rd b/man/svygof.Rd index 0225f25..39da871 100644 --- a/man/svygof.Rd +++ b/man/svygof.Rd @@ -20,8 +20,10 @@ that was used to fit the model.} fitted probabilities. Defaults to 10 (deciles).} } \value{ -A `data.frame` containing the F-statistic, the numerator (df1) and -denominator (df2) degrees of freedom, and the p-value for the test. +A one-row `data.frame` with columns `F_statistic`, `df1` (numerator df), +`df2` (denominator df) and `p_value` for the Archer-Lemeshow design-based +goodness-of-fit test. A non-significant p-value (e.g. p > 0.05) is consistent +with adequate fit; it does not prove the model is correctly specified. } \description{ Performs an Archer-Lemeshow goodness-of-fit (GOF) test for logistic @@ -36,41 +38,21 @@ Wald test. A non-significant p-value (e.g., p > 0.05) suggests no evidence of a poor fit. } \examples{ -\dontrun{ -# Ensure required packages are loaded -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE)) { +\donttest{ +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - # 1. Prepare Data - data(NHANESraw, package = "NHANES") - nhanes_data <- NHANESraw \%>\% - dplyr::filter(Age >= 20) \%>\% - dplyr::mutate(ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese"))) \%>\% - dplyr::filter(complete.cases(ObeseStatus, Age, Gender, Race1, - WTMEC2YR, SDMVPSU, SDMVSTRA)) +# The Archer-Lemeshow test uses a standard (non-replicate) survey design. +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # 2. Create a replicate design object - std_design <- survey::svydesign( - ids = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_data - ) - rep_design <- survey::as.svrepdesign(std_design) +fit <- survey::svyglm(htn01 ~ age + sex + smoking, + design = design, family = quasibinomial()) - # 3. Fit a survey logistic regression model using the replicate design - fit_obesity_rep <- survey::svyglm( - ObeseStatus ~ Age + Gender + Race1, - design = rep_design, - family = quasibinomial() - ) - - # 4. Calculate the design-correct AUC - auc_results <- svyAUC(fit_obesity_rep, rep_design) - print(auc_results) -} +# A non-significant p-value is consistent with adequate fit (it does not +# prove the model is correctly specified). +svygof(fit, design, G = 10) } } diff --git a/man/svykmplot.Rd b/man/svykmplot.Rd index 3016b3d..b1dd727 100644 --- a/man/svykmplot.Rd +++ b/man/svykmplot.Rd @@ -13,6 +13,7 @@ svykmplot( risk_table_title = "Number at Risk", palette = NULL, show_pval = TRUE, + se = TRUE, show_censor_marks = TRUE, base_font_size = 11, table_font_size = 3.5 @@ -40,6 +41,11 @@ use for the strata.} \item{show_pval}{Logical. If `TRUE`, calculates and displays the survey-weighted log-rank test p-value.} +\item{se}{Logical. If `TRUE` (default), computes survey-weighted confidence +bands via `survey::svykm(se = TRUE)` and draws them as shaded ribbons. The +standard-error calculation can be slow on large designs; set `se = FALSE` +to draw the survival curves only (much faster, no confidence bands).} + \item{show_censor_marks}{Logical. If `TRUE`, displays censor markings (`+`) on the survival curves.} @@ -62,53 +68,30 @@ with an attached "Number at Risk" table and censor markings, using `ggplot2` and `patchwork`. } \examples{ -\dontrun{ -# This example uses the 'nhanes_mortality' dataset, - -if (requireNamespace("survey", quietly = TRUE)) { - - # 1. Load the data - data(nhanes_mortality) - - # 2. Create the main survey design object - analytic_design <- survey::svydesign( - strata = ~strata, - id = ~psu, - weights = ~survey_weight, - data = nhanes_mortality, - nest = TRUE - ) - - # 3. Create a subsetted design for females - design_female <- subset(analytic_design, sex == "Female") - - # 4. Define the formula - km_formula <- Surv(stime, status) ~ caff +\donttest{ +library(survival) # for Surv() in the formula below +data(nhanes_mortality, package = "svyTable1") - # 5. Define a 4-color palette - distinct_palette <- c("#377EB8", "#FF7F00", "#4DAF4A", "#E41A1C") - - # 6. Run the function - km_results_female <- svykmplot( - formula = km_formula, - design = design_female, - legend_title = "Caffeine Consumption", - time_unit = "days", # Use "days" so divisor is 1 - time_breaks = seq(0, 240, by = 60), - palette = distinct_palette, - show_pval = TRUE, - show_censor_marks = TRUE - ) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # 7. Display the plot (and correct the x-axis label) - if (requireNamespace("ggplot2", quietly = TRUE)) { - km_results_female$plot + - ggplot2::labs(x = "Follow-up Time (Months)") - } +# Survey-weighted Kaplan-Meier curves by caffeine consumption (women only). +design_female <- subset(design, sex == "Female") - # 8. Print the at-risk table - print(km_results_female$table) +# se = FALSE draws curves only (fast); use se = TRUE to add confidence bands. +km <- svykmplot( + formula = Surv(stime, status) ~ caff, + design = design_female, + legend_title = "Caffeine consumption", + time_unit = "days", + time_breaks = seq(0, 240, by = 60), + show_pval = TRUE, + se = FALSE +) -} +km$plot # a ggplot of the weighted survival curves +print(km$table) # the number-at-risk table } } diff --git a/man/svypooled.Rd b/man/svypooled.Rd index 93ccb4f..1b8ddbc 100644 --- a/man/svypooled.Rd +++ b/man/svypooled.Rd @@ -34,8 +34,12 @@ returns a table showing only the main exposure and lists covariates in a footnote. If `FALSE`, it returns a full table with all model terms.} } \value{ -An HTML table object of class `kableExtra`. This object can be - printed directly in R Markdown documents or saved. +A `kable` object (class `c("kableExtra", "knitr_kable")`): a character + vector of HTML with attributes, ready to print in an R Markdown report. In + "fallacy-safe" mode it contains one row per level of `main_exposure` (with the + adjustment variables named in a footnote); otherwise it contains all model + terms grouped by variable. Each cell shows the effect measure and its 95\% + confidence interval pooled across imputations using Rubin's rules. } \description{ This function takes a pooled model object from the `mice` package and creates a @@ -58,62 +62,31 @@ formats the results into a clean HTML table using `knitr::kable()` and less than 0.001 shown as "<0.001". } \examples{ -\dontrun{ -# Load required packages for the example -library(mice) -library(dplyr) -library(survey) -library(NHANES) +\donttest{ +data(nhanes_mortality, package = "svyTable1") +dat <- nhanes_mortality +dat$htn01 <- as.numeric(dat$htn == "Yes") +dat <- dat[, c("htn01", "sex", "age", "smoking", "psu", "strata", "survey_weight")] -# --- 1. Data Preparation --- -# We'll prepare a clean analytic dataset from the raw NHANES data. -# Note: We convert the outcome 'Obese' to a numeric 0/1 variable for svyglm. -data(NHANESraw, package = "NHANES") -nhanes_analytic <- NHANESraw \%>\% - dplyr::filter(Age >= 20 & WTMEC2YR > 0) \%>\% - mutate( - Obese_numeric = ifelse(BMI >= 30, 1, 0), - AgeCat = cut(Age, breaks = c(19, 39, 59, 80), labels = c("20-39", "40-59", "60-80")), - Smoke100 = factor(Smoke100, levels = c("No", "Yes")) - ) \%>\% - select(Obese_numeric, AgeCat, Smoke100, Education, SDMVPSU, SDMVSTRA, WTMEC2YR) +# Induce missingness, then multiply impute (small m for a fast example only; +# for a real analysis use m >= 20 and check convergence). +set.seed(1) +dat$age[sample(nrow(dat), 200)] <- NA +imp <- mice::mice(dat, m = 2, maxit = 2, printFlag = FALSE, seed = 1) -# --- 2. Perform Imputation and Pooled Analysis --- -# Set survey option to handle lonely PSUs, a common issue with NHANES data. -options(survey.lonely.psu = "adjust") +# Fit a survey-weighted model in each imputed dataset, then pool (Rubin's rules). +fits <- with(imp, survey::svyglm( + htn01 ~ sex + age + smoking, + design = survey::svydesign(id = ~psu, strata = ~strata, + weights = ~survey_weight, nest = TRUE, + data = data.frame(mget(ls()))), + family = quasibinomial() +)) +pooled <- mice::pool(fits) -# Impute the analytic dataset -imputed_data <- mice(nhanes_analytic, m = 2, maxit = 2, seed = 123, printFlag = FALSE) - -# Use with() to run a survey-weighted GLM on each imputed dataset -fit_imputed <- with(imputed_data, - svyglm(Obese_numeric ~ Smoke100 + AgeCat + Education, - design = svydesign(id = ~SDMVPSU, strata = ~SDMVSTRA, - weights = ~WTMEC2YR, nest = TRUE, data = .data), - family = quasibinomial()) - ) - -# Pool the results from all models into a single 'mipo' object -pooled_results <- pool(fit_imputed) - -# --- 3. Generate Tables with svypooled --- -# Example A: Create a "fallacy-safe" table (default) -svypooled( - pooled_model = pooled_results, - main_exposure = "Smoke100", - adj_var_names = c("AgeCat", "Education"), - measure = "OR", - title = "Adjusted Odds of Obesity (Fallacy-Safe)" -) - -# Example B: Create a full table with all variables -svypooled( - pooled_model = pooled_results, - main_exposure = "Smoke100", - adj_var_names = c("AgeCat", "Education"), - measure = "OR", - title = "Adjusted Odds of Obesity (Full Table for Appendix)", - fallacy_safe = FALSE -) +# Fallacy-safe table: shows only the exposure; covariates go in a footnote. +svypooled(pooled, main_exposure = "sex", + adj_var_names = c("age", "smoking"), measure = "OR", + title = "Adjusted odds of hypertension") } } diff --git a/man/svytable1.Rd b/man/svytable1.Rd index 89107d9..f0754cb 100644 --- a/man/svytable1.Rd +++ b/man/svytable1.Rd @@ -51,8 +51,18 @@ reliability metrics. This includes the RSE and specific TRUE/FALSE columns (e.g., `fail_n_30`, `fail_eff_n_30`) for each suppression rule. Defaults to FALSE.} } \value{ -If `return_metrics` is FALSE (default), returns a data.frame. If TRUE, - returns a list with two elements: `formatted_table` and `reliability_metrics`. +If `return_metrics` is FALSE (default), a `data.frame` (the formatted + Table 1). Its first two columns are `Variable` and `Level`; the remaining + columns are `Overall` plus one column per level of `strata_var` (including a + `Missing` column when the stratifier has NAs). Cells hold unweighted counts + with weighted percentages for categorical variables and "mean (SD)" for + numeric variables; suppressed cells (when `reliability_checks = TRUE`) are + shown as `"*"`. If `return_metrics` is TRUE, a list with two elements: + `formatted_table` (the data.frame above) and `reliability_metrics` (a + data.frame with one row per evaluated cell giving `n`, `df`, `deff`, + `effective_n`, the confidence-interval bounds, `rse`, the logical + suppression flags `fail_n_30`/`fail_eff_n_30`/`fail_df_8`/`fail_ciw_30`/ + `fail_rciw_130`/`fail_rse_30`, and overall `suppressed`). } \description{ Generates a "Table 1" of descriptive statistics for complex survey data, diff --git a/tests/testthat/helper-svy.R b/tests/testthat/helper-svy.R new file mode 100644 index 0000000..9c507a8 --- /dev/null +++ b/tests/testthat/helper-svy.R @@ -0,0 +1,38 @@ +# ---- Shared fixtures for svyTable1 tests ---------------------------------- + +# Deterministic mock model exposing coef()/vcov(), used for the pure-math +# (Tier 1) additive-interaction tests so the expected numbers do not depend on +# the survey package's variance machinery. +coef.svymock <- function(object, ...) object$coefficients +vcov.svymock <- function(object, ...) object$vcov +.S3method("coef", "svymock", coef.svymock) +.S3method("vcov", "svymock", vcov.svymock) + +make_mock <- function(beta, V) { + if (is.null(dimnames(V))) dimnames(V) <- list(names(beta), names(beta)) + structure(list(coefficients = beta, vcov = V), class = "svymock") +} + +# Analytic data + standard design on the BUNDLED dataset (no Suggests needed). +st1_data <- function() { + e <- new.env() + utils::data("nhanes_mortality", package = "svyTable1", envir = e) + d <- e$nhanes_mortality + d$htn01 <- as.numeric(d$htn == "Yes") + d +} + +st1_design <- function(data = st1_data()) { + survey::svydesign(id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = data) +} + +st1_fit <- function(design = st1_design()) { + survey::svyglm(htn01 ~ age + sex + smoking, design = design, + family = stats::quasibinomial()) +} + +st1_int_fit <- function(design = st1_design()) { + survey::svyglm(htn01 ~ sex * insulin + age, design = design, + family = stats::quasibinomial()) +} diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R new file mode 100644 index 0000000..0716a57 --- /dev/null +++ b/tests/testthat/setup.R @@ -0,0 +1,3 @@ +# Survey SEs for NHANES-style designs require an explicit lonely-PSU policy. +# Tests run in their own process, so setting it here does not leak to the user. +options(survey.lonely.psu = "adjust") diff --git a/tests/testthat/test-addint.R b/tests/testthat/test-addint.R new file mode 100644 index 0000000..d7f57a5 --- /dev/null +++ b/tests/testthat/test-addint.R @@ -0,0 +1,78 @@ +# Tier 1: pure-math additive interaction on FIXED coef/vcov. These numbers are +# fully deterministic (no survey variance machinery), so they are pinned tightly. + +test_that("addint computes RERI/AP/S exactly (interaction parameterisation)", { + beta <- c(A = 0.4, B = 0.3, `A:B` = 0.5) + V <- diag(c(0.02, 0.03, 0.05)) + dimnames(V) <- list(names(beta), names(beta)) + m <- make_mock(beta, V) + + res <- addint(m, type = "interaction", + coef_names = list(exp1_coef = "A", exp2_coef = "B", + inter_coef = "A:B"), + measures = "all") + + hr10 <- exp(beta[["A"]]); hr01 <- exp(beta[["B"]]); hr11 <- exp(sum(beta)) + reri <- hr11 - hr10 - hr01 + 1 + ap <- reri / hr11 + s <- (hr11 - 1) / ((hr10 - 1) + (hr01 - 1)) + + expect_equal(unname(res$RERI[["RERI_Estimate"]]), reri, tolerance = 1e-8) + expect_equal(unname(res$AP[["AP_Estimate"]]), ap, tolerance = 1e-8) + expect_equal(unname(res$S[["S_Estimate"]]), s, tolerance = 1e-8) + + # CI is a symmetric Wald interval around the point estimate on its scale. + z <- stats::qnorm(0.975) + expect_equal(unname(res$RERI[["RERI_CI95_low"]]), + reri - z * unname(res$RERI[["RERI_SE"]]), tolerance = 1e-8) + expect_equal(unname(res$RERI[["RERI_CI95_upp"]]), + reri + z * unname(res$RERI[["RERI_SE"]]), tolerance = 1e-8) + + # Delta-method SEs are deterministic: pin them as a regression guard. + expect_equal(unname(res$RERI[["RERI_SE"]]), 0.8570119, tolerance = 1e-5) + expect_equal(unname(res$AP[["AP_SE"]]), 0.1275353, tolerance = 1e-5) +}) + +test_that("addint joint and interaction parameterisations give the same point estimate", { + beta <- c(A = 0.4, B = 0.3, `A:B` = 0.5) + b_joint <- c(a10 = beta[["A"]], a01 = beta[["B"]], a11 = sum(beta)) + Vj <- diag(c(0.02, 0.03, 0.08)) + dimnames(Vj) <- list(names(b_joint), names(b_joint)) + mj <- make_mock(b_joint, Vj) + + rj <- addint(mj, type = "joint", + coef_names = list(exp1_level = "a10", exp2_level = "a01", + both_levels = "a11"), + measures = "RERI") + + hr10 <- exp(0.4); hr01 <- exp(0.3); hr11 <- exp(1.2) + expect_equal(unname(rj$RERI[["RERI_Estimate"]]), + hr11 - hr10 - hr01 + 1, tolerance = 1e-8) +}) + +test_that("addint warns and returns NULL when a required coefficient is missing", { + beta <- c(A = 0.4, B = 0.3) + V <- diag(2) + dimnames(V) <- list(names(beta), names(beta)) + m <- make_mock(beta, V) + expect_warning( + res <- addint(m, type = "interaction", + coef_names = list(exp1_coef = "A", exp2_coef = "B", + inter_coef = "A:B")), + "not found" + ) + expect_null(res) +}) + +test_that("addint validates the measures argument", { + beta <- c(A = 0.4, B = 0.3, `A:B` = 0.5) + V <- diag(3) + dimnames(V) <- list(names(beta), names(beta)) + m <- make_mock(beta, V) + expect_error( + addint(m, type = "interaction", + coef_names = list(exp1_coef = "A", exp2_coef = "B", inter_coef = "A:B"), + measures = "bogus"), + "Invalid measure" + ) +}) diff --git a/tests/testthat/test-diagnostics.R b/tests/testthat/test-diagnostics.R new file mode 100644 index 0000000..2dbcc5a --- /dev/null +++ b/tests/testthat/test-diagnostics.R @@ -0,0 +1,63 @@ +# svydiag -------------------------------------------------------------------- + +test_that("svydiag returns a tibble with the documented columns", { + out <- svydiag(st1_fit()) + expect_s3_class(out, "tbl_df") + expect_true(all(c("Term", "Estimate", "SE", "p.value", "is_significant", + "CI_Lower", "CI_Upper", "CI_Width", "RSE_percent", + "is_rse_high") %in% names(out))) + expect_type(out$is_significant, "logical") + expect_true(all(out$CI_Width >= 0)) +}) + +test_that("svydiag honours custom thresholds", { + out <- svydiag(st1_fit(), p_threshold = 1) + expect_true(all(out$is_significant)) +}) + +# svygof --------------------------------------------------------------------- + +test_that("svygof returns a one-row F-test on a standard design", { + des <- st1_design() + g <- svygof(st1_fit(des), des, G = 10) + expect_s3_class(g, "data.frame") + expect_equal(nrow(g), 1) + expect_named(g, c("F_statistic", "df1", "df2", "p_value")) + expect_gte(g$p_value, 0) + expect_lte(g$p_value, 1) +}) + +test_that("svygof works on a replicate design (regression: it used to error)", { + des <- st1_design() + rep_des <- survey::as.svrepdesign(des) + fit <- survey::svyglm(htn01 ~ age + sex + smoking, design = rep_des, + family = stats::quasibinomial()) + g <- svygof(fit, rep_des, G = 10) + expect_s3_class(g, "data.frame") + expect_gte(g$p_value, 0) + expect_lte(g$p_value, 1) +}) + +test_that("svygof rejects a non-svyglm fit", { + expect_error(svygof(stats::lm(mpg ~ wt, data = mtcars), st1_design()), "svyglm") +}) + +# svyAUC --------------------------------------------------------------------- + +test_that("svyAUC returns a CI bounded within [0, 1]", { + des <- st1_design() + rep_des <- survey::as.svrepdesign(des) + fit <- survey::svyglm(htn01 ~ age + sex + smoking, design = rep_des, + family = stats::quasibinomial()) + a <- svyAUC(fit, rep_des) + expect_s3_class(a, "data.frame") + expect_named(a, c("AUC", "SE", "CI_Lower", "CI_Upper")) + expect_gte(a$CI_Lower, 0) + expect_lte(a$CI_Upper, 1) + expect_true(a$AUC >= a$CI_Lower && a$AUC <= a$CI_Upper) +}) + +test_that("svyAUC requires a replicate-weights design", { + des <- st1_design() + expect_error(svyAUC(st1_fit(des), des), "replicate") +}) diff --git a/tests/testthat/test-interaction.R b/tests/testthat/test-interaction.R new file mode 100644 index 0000000..7a9fad4 --- /dev/null +++ b/tests/testthat/test-interaction.R @@ -0,0 +1,50 @@ +# Higher-level interaction reporting on the bundled survey design. + +test_that("jointeffects returns one row per factor-level combination", { + je <- jointeffects(st1_int_fit(), "sex", "insulin") + expect_s3_class(je, "tbl_df") + expect_true(all(c("Level1", "Level2", "Estimate", "SE", "CI.low", "CI.upp") %in% + names(je))) + expect_equal(nrow(je), 4) # 2 x 2 factor combinations +}) + +test_that("addintlist returns the three additive-interaction measures", { + al <- addintlist(model = st1_int_fit(), factor1_name = "sex", + factor2_name = "insulin", measures = "all") + expect_s3_class(al, "tbl_df") + expect_true(all(c("Factor1", "Level1", "Factor2", "Level2", "Measure", + "Estimate", "SE", "CI_low", "CI_upp") %in% names(al))) + expect_setequal(unique(al$Measure), c("RERI", "AP", "S")) +}) + +test_that("reportint(output = 'list') returns the documented panels", { + fit <- st1_int_fit() + suppressMessages(invisible(utils::capture.output( + rl <- reportint(fit, factor1_name = "sex", factor2_name = "insulin", + output = "list") + ))) + expect_type(rl, "list") + expect_true(all(c("joint_effects", "stratum_specific_effects", + "additive_interaction", "multiplicative_scale", + "effect_modification_report", "Interaction_report") %in% + names(rl))) + expect_s3_class(rl$joint_effects, "data.frame") +}) + +test_that("reportint auto-detects a single interaction term", { + fit <- st1_int_fit() + suppressMessages(invisible(utils::capture.output( + rl <- reportint(fit, output = "list") + ))) + expect_type(rl, "list") +}) + +test_that("reportint errors when no interaction term is present", { + des <- st1_design() + fit <- survey::svyglm(htn01 ~ sex + insulin + age, design = des, + family = stats::quasibinomial()) + expect_error( + suppressMessages(reportint(fit, output = "list")), + "interaction" + ) +}) diff --git a/tests/testthat/test-smoke.R b/tests/testthat/test-smoke.R deleted file mode 100644 index 6c46793..0000000 --- a/tests/testthat/test-smoke.R +++ /dev/null @@ -1,9 +0,0 @@ -test_that("svytable1 runs without error", { - # A very simple check to ensure the function is available and runs - expect_true(is.function(svytable1)) -}) - -test_that("svypooled runs without error", { - # A simple check for the other main function - expect_true(is.function(svypooled)) -}) diff --git a/tests/testthat/test-survival.R b/tests/testthat/test-survival.R new file mode 100644 index 0000000..cb7af10 --- /dev/null +++ b/tests/testthat/test-survival.R @@ -0,0 +1,37 @@ +# svykmplot ------------------------------------------------------------------ + +test_that("svykmplot returns a plot and an at-risk table", { + des <- st1_design() + des_f <- subset(des, sex == "Female") + # se = FALSE keeps this fast; the confidence-band path (se = TRUE) is the + # long-standing default and is exercised separately below on a small subset. + km <- svykmplot(survival::Surv(stime, status) ~ caff, design = des_f, + time_unit = "days", time_breaks = seq(0, 240, by = 60), + show_pval = TRUE, se = FALSE) + expect_type(km, "list") + expect_true(all(c("plot", "table") %in% names(km))) + expect_s3_class(km$plot, "ggplot") +}) + +test_that("svykmplot draws confidence bands when se = TRUE", { + des <- st1_design() + # A small subgroup keeps the (slow) svykm SE computation quick. + des_small <- subset(des, sex == "Female" & race == "Mexican American") + km <- svykmplot(survival::Surv(stime, status) ~ caff, design = des_small, + time_unit = "days", time_breaks = seq(0, 240, by = 60), + show_pval = FALSE, se = TRUE) + expect_s3_class(km$plot, "ggplot") +}) + +# svycoxph_CE ---------------------------------------------------------------- + +test_that("svycoxph_CE returns a ggplot constant-effect diagnostic", { + des <- st1_design() + suppressMessages(invisible(utils::capture.output( + p <- svycoxph_CE(formula_rhs = "sex + age", design = des, + var_to_test = "sexFemale", time_var = "stime", + status_var = "status", n_intervals = 3, + print_main_model = FALSE, print_split_summary = FALSE) + ))) + expect_s3_class(p, "ggplot") +}) diff --git a/tests/testthat/test-svypooled.R b/tests/testthat/test-svypooled.R new file mode 100644 index 0000000..1c17dcd --- /dev/null +++ b/tests/testthat/test-svypooled.R @@ -0,0 +1,43 @@ +# svypooled renders a publication table from a mice::pool() (mipo) object. + +make_pooled <- function() { + d <- st1_data() + da <- d[, c("htn01", "sex", "age", "smoking", "psu", "strata", "survey_weight")] + set.seed(1) + da$age[sample(nrow(da), 200)] <- NA + imp <- mice::mice(da, m = 2, maxit = 2, printFlag = FALSE, seed = 1) + fit_imp <- with(imp, survey::svyglm( + htn01 ~ sex + age + smoking, + design = survey::svydesign(id = ~psu, strata = ~strata, + weights = ~survey_weight, nest = TRUE, + data = data.frame(mget(ls()))), + family = stats::quasibinomial() + )) + mice::pool(fit_imp) +} + +test_that("svypooled returns a kableExtra/knitr_kable object (fallacy-safe)", { + skip_if_not_installed("mice") + tab <- svypooled(make_pooled(), main_exposure = "sex", + adj_var_names = c("age", "smoking"), measure = "OR", + title = "Test") + expect_s3_class(tab, "knitr_kable") + expect_s3_class(tab, "kableExtra") +}) + +test_that("svypooled full table also renders", { + skip_if_not_installed("mice") + tab <- svypooled(make_pooled(), main_exposure = "sex", + adj_var_names = c("age", "smoking"), measure = "OR", + fallacy_safe = FALSE) + expect_s3_class(tab, "knitr_kable") +}) + +test_that("svypooled errors when the main exposure is absent", { + skip_if_not_installed("mice") + expect_error( + svypooled(make_pooled(), main_exposure = "not_there", + adj_var_names = c("age", "smoking")), + "not found" + ) +}) diff --git a/tests/testthat/test-svytable1.R b/tests/testthat/test-svytable1.R new file mode 100644 index 0000000..e700df7 --- /dev/null +++ b/tests/testthat/test-svytable1.R @@ -0,0 +1,73 @@ +test_that("svytable1 returns a data.frame with the expected structure", { + des <- st1_design() + tab <- svytable1(des, strata_var = "htn", table_vars = c("age", "sex", "smoking")) + expect_s3_class(tab, "data.frame") + expect_true(all(c("Variable", "Level", "Overall", "No", "Yes") %in% names(tab))) + expect_equal(tab$Variable[1], "n") +}) + +test_that("svytable1 mixed-mode N row reports unweighted sample sizes", { + d <- st1_data() + des <- st1_design(d) + tab <- svytable1(des, "htn", "age") + # Default "mixed" mode: the Overall n is the real (unweighted) row count. + expect_equal(tab$Overall[1], format(nrow(d), big.mark = ",")) +}) + +test_that("svytable1 errors on variables not in the design", { + des <- st1_design() + expect_error(svytable1(des, "htn", c("age", "not_a_var")), "not found") +}) + +test_that("svytable1 errors loudly on a 0-row design", { + des <- st1_design() + empty <- subset(des, age > 1e6) + expect_error(svytable1(empty, "htn", "age"), "0 rows") +}) + +test_that("svytable1 weighted and unweighted modes run", { + des <- st1_design() + expect_s3_class(svytable1(des, "htn", "age", mode = "weighted"), "data.frame") + expect_s3_class(svytable1(des, "htn", "age", mode = "unweighted"), "data.frame") +}) + +test_that("svytable1 reliability metrics expose the documented NCHS flags", { + des <- st1_design() + res <- svytable1(des, "htn", c("age", "sex", "smoking"), + reliability_checks = TRUE, return_metrics = TRUE) + expect_type(res, "list") + expect_named(res, c("formatted_table", "reliability_metrics")) + expect_true(all(c("fail_n_30", "fail_eff_n_30", "fail_df_8", "fail_ciw_30", + "fail_rciw_130", "fail_rse_30", "suppressed") %in% + names(res$reliability_metrics))) + expect_type(res$reliability_metrics$suppressed, "logical") +}) + +test_that("svytable1 suppresses a deliberately tiny, unreliable cell", { + # Build a tiny stratum so at least one categorical cell fails NCHS rules. + d <- st1_data() + des <- st1_design(d) + res <- svytable1(des, "htn", c("race"), + reliability_checks = TRUE, return_metrics = TRUE) + expect_true(any(res$reliability_metrics$suppressed)) + # A suppressed cell renders as "*" in the formatted table. + expect_true(any(vapply(res$formatted_table, function(col) any(col == "*"), logical(1)))) +}) + +test_that("svytable1 folds NA in the stratifier into a 'Missing' column", { + d <- st1_data() + set.seed(11) + d$htn[sample(nrow(d), 50)] <- NA + des <- st1_design(d) + tab <- svytable1(des, "htn", "age") + expect_true("Missing" %in% names(tab)) +}) + +test_that("svytable1 reports a Missing row for a numeric variable with NAs", { + d <- st1_data() + set.seed(12) + d$carbohyd[sample(nrow(d), 100)] <- NA + des <- st1_design(d) + tab <- svytable1(des, "htn", "carbohyd") + expect_true(any(grepl("Missing", tab$Level))) +})