diff --git a/DESCRIPTION b/DESCRIPTION index 1796b0d..4d2b5ba 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -57,7 +57,6 @@ Imports: WeightedROC, broom Suggests: - NHANES, mitools, testthat (>= 3.0.0) VignetteBuilder: knitr diff --git a/R/addint.R b/R/addint.R index 255e0f7..a57e164 100644 --- a/R/addint.R +++ b/R/addint.R @@ -22,6 +22,15 @@ #' #' Confidence intervals for S are calculated on the log scale and then exponentiated. #' +#' \strong{Scale caveat:} these measures are defined on the risk-ratio scale. +#' When the model supplies odds ratios (logistic) or hazard ratios (Cox), RERI, +#' AP, and S only approximate their risk-ratio counterparts if the outcome is +#' rare; for a common outcome the odds ratio overstates the risk ratio and the +#' measures can be biased. To target risk ratios directly, fit a log-binomial or +#' Poisson working model (see VanderWeele and Knol (2014) +#' ). The delta-method Wald intervals here are +#' symmetric and approximate. +#' #' @param model A fitted model object (e.g., `svycoxph`, `svyglm`). #' @param type Character string: `"joint"` if using a combined categorical variable, #' `"interaction"` if using main effects and a product term. diff --git a/R/addintlist.R b/R/addintlist.R index 80c2743..cdf5f28 100644 --- a/R/addintlist.R +++ b/R/addintlist.R @@ -40,9 +40,21 @@ #' \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. +#' \item `Scale`: The ratio scale the measures were computed on, inferred +#' from the model: "OR (logistic)", "HR (Cox)", "RR (log-binomial)", +#' "RR (Poisson)", or "ratio". #' } #' Returns an empty tibble if errors occur. #' +#' @details +#' RERI, AP, and the synergy index S are additive-interaction measures defined on +#' the \emph{risk-ratio} scale. When the model reports odds ratios (logistic) or +#' hazard ratios (Cox), these measures only approximate their risk-ratio +#' counterparts when the outcome is rare; for a common outcome the odds ratio +#' overstates the risk ratio and RERI/AP/S can be biased. The `Scale` column makes +#' the operating scale explicit. To estimate risk-ratio-based interaction +#' directly for a common outcome, fit a log-binomial or Poisson working model. +#' #' @seealso \code{\link{addint}} #' #' @importFrom stats coef @@ -55,65 +67,28 @@ #' #' @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() +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' adult_design_binary <- svydesign( -#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR, -#' nest = TRUE, data = nhanes_adults_processed -#' ) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # --- 2. Fit Interaction Term Model --- -#' interaction_model_logit <- svyglm( -#' Hypertension_130 ~ Race1 * ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) +#' # Saturated interaction model of sex and insulin use. +#' interaction_model <- survey::svyglm( +#' htn01 ~ sex * insulin + age, +#' design = design, 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.") -#' } +#' # RERI, AP and S for every factor-level combination. The Scale column flags +#' # that these are OR-based (see Details for the rare-outcome caveat). +#' addintlist( +#' model = interaction_model, +#' factor1_name = "sex", +#' factor2_name = "insulin", +#' measures = "all" +#' ) #' } addintlist <- function(model, factor1_name, @@ -307,10 +282,29 @@ addintlist <- function(model, } # --- Final structure --- + # The Scale column makes explicit which ratio measure RERI/AP/S were built + # from (OR for logistic, HR for Cox, RR for log-binomial/Poisson), since these + # additive-interaction measures are only RR-scale-correct when the outcome is + # rare (see Details). final_table <- final_table %>% - # Use .data pronoun + dplyr::mutate(Scale = .effect_measure_scale(model)) %>% dplyr::relocate("Factor1", "Level1", "Factor2", "Level2", - "Measure", "Estimate", "SE", "CI_low", "CI_upp") + "Measure", "Estimate", "SE", "CI_low", "CI_upp", "Scale") return(final_table) } + +# Internal: infer the ratio scale that RERI/AP/S were computed on, from the +# fitted model's class/family. Not exported. +.effect_measure_scale <- function(model) { + cl <- class(model)[1] + if (cl %in% c("svycoxph", "coxph")) return("HR (Cox)") + fam <- tryCatch(stats::family(model)$family, error = function(e) NA_character_) + link <- tryCatch(stats::family(model)$link, error = function(e) NA_character_) + if (!is.na(fam)) { + if (grepl("binomial", fam) && identical(link, "logit")) return("OR (logistic)") + if (grepl("binomial", fam) && identical(link, "log")) return("RR (log-binomial)") + if (grepl("poisson", fam)) return("RR (Poisson)") + } + "ratio" +} diff --git a/R/jointeffects.R b/R/jointeffects.R index 3060f9e..2a9931a 100644 --- a/R/jointeffects.R +++ b/R/jointeffects.R @@ -47,70 +47,26 @@ #' #' @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 ratio SE +#' 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 +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # --- 1. Data Preparation (NHANES Example) --- -#' data(NHANESraw) +#' interaction_model <- survey::svyglm( +#' htn01 ~ sex * insulin + age, +#' design = design, family = quasibinomial() +#' ) #' -#' 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 the Interaction Model --- -#' interaction_model_logit <- svyglm( -#' Hypertension_130 ~ Race1 * ObeseStatus + Age, -#' design = adult_design_binary, family = quasibinomial() -#' ) -#' -#' # --- 3. Run the function --- -#' -#' # Example 1: Output on Ratio scale (default) -#' joint_effects_ratio <- jointeffects( -#' interaction_model = interaction_model_logit, -#' factor1_name = "Race1", -#' factor2_name = "ObeseStatus", -#' scale = "ratio" -#' ) -#' print("--- Joint Effects (Ratio Scale) ---") -#' print(joint_effects_ratio, n = 50) -#' -#' # Example 2: Output on Log scale -#' joint_effects_log <- jointeffects( -#' interaction_model = interaction_model_logit, -#' factor1_name = "Race1", -#' factor2_name = "ObeseStatus", -#' scale = "log" -#' ) -#' print("--- Joint Effects (Log Scale) ---") -#' print(joint_effects_log, n = 50) -#' } +#' # Joint effects of sex x insulin on the ratio (OR) scale. +#' jointeffects( +#' interaction_model = interaction_model, +#' factor1_name = "sex", +#' factor2_name = "insulin", +#' scale = "ratio" +#' ) #' } jointeffects <- function(interaction_model, factor1_name, diff --git a/R/plotint.R b/R/plotint.R index 1865b9b..73d6faa 100644 --- a/R/plotint.R +++ b/R/plotint.R @@ -74,95 +74,8 @@ #' legend_title = "log(Total Calories)" #' ) #' -#' # You can now access the plot data if you want it -#' # head(plot_data$plot_data) -#' -#' -#' # ================================================================= -#' # Example 2: Continuous x Continuous (NHANESraw) -#' # ================================================================= -#' library(NHANES) -#' library(dplyr) -#' library(survey) -#' library(emmeans) -#' library(ggplot2) -#' -#' # --- 3. DATA PREPARATION (NHANESraw) --- -#' data("NHANESraw", package = "NHANES") -#' -#' # Select relevant variables and filter for adults with complete cases -#' nhanes_clean <- NHANESraw %>% -#' filter(Age >= 20) %>% -#' mutate( -#' # Ensure factors are set correctly -#' Gender = factor(Gender), -#' Race1 = factor(Race1) -#' ) %>% -#' # Keep only the variables we need for this model -#' select( -#' BPSysAve, Age, BMI, Gender, Race1, -#' SDMVPSU, SDMVSTRA, WTMEC2YR -#' ) %>% -#' na.omit() # Listwise deletion for simplicity -#' -#' # --- 4. DESIGN & MODEL FITTING --- -#' # Create the survey design object -#' design_nhanes <- svydesign( -#' id = ~SDMVPSU, -#' strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, -#' nest = TRUE, -#' data = nhanes_clean -#' ) -#' -#' # Fit the model: Interaction of Age * BMI on Geometric Mean of BP -#' # We use gaussian(link="log") to model Geometric Means -#' model_bp_interaction <- svyglm( -#' BPSysAve ~ Age * BMI + Gender + Race1, -#' design = design_nhanes, -#' family = gaussian(link = "log") -#' ) -#' -#' # --- 5. PLOT THE INTERACTION --- -#' -#' # --- Plot 1: Standard Plot (using quantiles and CIs) --- -#' # This shows the effect of Age (x-axis) at 10th/50th/90th percentiles of BMI -#' print( -#' plotint( -#' model = model_bp_interaction, -#' effect = "Age", -#' moderator = "BMI", -#' data = nhanes_clean, -#' mod_scale = "quantile", # Use 10th/50th/90th percentiles for BMI -#' type = "response", # Plot the Geometric Mean (back-transform from log) -#' show_ci = TRUE, -#' bw = FALSE, -#' pval_position = "bottom.right", -#' xlab = "Age (Years)", -#' ylab = "Predicted Geometric Mean SBP (mmHg)", -#' legend_title = "BMI (Percentiles)" -#' ) -#' ) -#' -#' # --- Plot 2: Cleaner Plot (using SD and CIs) --- -#' # This shows the effect of Age at Mean +/- 1 SD of BMI -#' # We turn off CIs for a cleaner look if they are too wide -#' print( -#' plotint( -#' model = model_bp_interaction, -#' effect = "Age", -#' moderator = "BMI", -#' data = nhanes_clean, -#' mod_scale = "sd", # Use Mean +/- 1 SD for BMI -#' type = "response", -#' show_ci = FALSE, # Turn off CI ribbons -#' bw = TRUE, # Use B&W linetypes -#' pval_position = "bottom.right", -#' xlab = "Age (Years)", -#' ylab = "Predicted Geometric Mean SBP (mmHg)", -#' legend_title = "BMI (Standard Deviations)" -#' ) -#' ) +#' # You can now access the underlying plot data: +#' head(plot_data$plot_data) #' } plotint <- function(model, effect, diff --git a/R/svydiag.R b/R/svydiag.R index 0ff0ae1..ca61bd3 100644 --- a/R/svydiag.R +++ b/R/svydiag.R @@ -48,54 +48,19 @@ #' @export #' #' @examples -#' # Ensure required packages are loaded -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE)) { +#' data(nhanes_mortality, package = "svyTable1") +#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") #' -#' # 1. Prepare Data using the NHANES example -#' data(NHANESraw, package = "NHANES") -#' nhanes_adults_with_na <- NHANESraw %>% -#' dplyr::filter(Age >= 20) %>% -#' dplyr::mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")), -#' Race1 = factor(Race1) -#' ) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' # Create a complete-case design object for the regression model -#' nhanes_complete <- nhanes_adults_with_na[complete.cases( -#' nhanes_adults_with_na[, c("ObeseStatus", "Age", "Race1")] -#' ), ] +#' fit <- survey::svyglm(htn01 ~ age + sex + race, +#' design = design, family = quasibinomial()) #' -#' adult_design_complete <- survey::svydesign( -#' id = ~SDMVPSU, -#' strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, -#' nest = TRUE, -#' data = nhanes_complete -#' ) -#' -#' # 2. Fit a survey-weighted logistic regression model -#' fit <- survey::svyglm( -#' ObeseStatus ~ Age + Race1, -#' design = adult_design_complete, -#' family = quasibinomial() -#' ) -#' -#' # 3. Get the reliability diagnostics table using the new function -#' diagnostics_table <- svydiag(fit) -#' -#' # Print the resulting table -#' print(diagnostics_table) -#' -#' # For a publication-ready table, pipe the result to kable() -#' if (requireNamespace("knitr", quietly = TRUE)) { -#' knitr::kable(diagnostics_table, -#' caption = "Reliability Diagnostics for NHANES Obesity Model", -#' digits = 3) -#' } -#' } +#' # Per-coefficient reliability diagnostics (SE, p-value, CI width, RSE). +#' svydiag(fit) svydiag <- function(fit, p_threshold = 0.05, rse_threshold = 30) { diff --git a/R/svytable1.R b/R/svytable1.R index 3576d1f..3aeccf0 100644 --- a/R/svytable1.R +++ b/R/svytable1.R @@ -54,52 +54,31 @@ #' @export #' #' @examples -#' # Ensure required packages are loaded -#' if (requireNamespace("survey", quietly = TRUE) && -#' requireNamespace("NHANES", quietly = TRUE) && -#' requireNamespace("dplyr", quietly = TRUE)) { +#' data(nhanes_mortality, package = "svyTable1") #' -#' # 1. Prepare Data using NHANES -#' data(NHANESraw, package = "NHANES") -#' nhanes_adults_with_na <- NHANESraw %>% -#' dplyr::filter(Age >= 20) %>% -#' dplyr::mutate( -#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), -#' levels = c("Not Obese", "Obese")) -#' ) +#' design <- survey::svydesign( +#' id = ~psu, strata = ~strata, weights = ~survey_weight, +#' nest = TRUE, data = nhanes_mortality +#' ) #' -#' adult_design_with_na <- survey::svydesign( -#' id = ~SDMVPSU, -#' strata = ~SDMVSTRA, -#' weights = ~WTMEC2YR, -#' nest = TRUE, -#' data = nhanes_adults_with_na -#' ) +#' # Basic survey-weighted Table 1, stratified by hypertension status. +#' svytable1( +#' design = design, +#' strata_var = "htn", +#' table_vars = c("age", "sex", "race", "smoking") +#' ) #' -#' # 2. Basic Example: Create a simple Table 1 -#' vars_to_summarize <- c("Age", "Gender", "Race1", "Education") -#' table1 <- svytable1( -#' design = adult_design_with_na, -#' strata_var = "ObeseStatus", -#' table_vars = vars_to_summarize -#' ) -#' print(table1) -#' -#' # 3. Advanced Example: Use reliability checks and get detailed metrics -#' results <- svytable1( -#' design = adult_design_with_na, -#' strata_var = "ObeseStatus", -#' table_vars = vars_to_summarize, -#' reliability_checks = TRUE, -#' return_metrics = TRUE -#' ) -#' -#' # View the table with unreliable estimates suppressed (*) -#' print(results$formatted_table) -#' -#' # View the detailed report card of reliability checks -#' print(results$reliability_metrics) -#' } +#' # With NCHS reliability checks plus a detailed metrics table. Estimates that +#' # fail the reliability rules are shown as "*". +#' results <- svytable1( +#' design = design, +#' strata_var = "htn", +#' table_vars = c("age", "sex", "race", "smoking"), +#' reliability_checks = TRUE, +#' return_metrics = TRUE +#' ) +#' results$formatted_table +#' head(results$reliability_metrics) svytable1 <- function(design, strata_var, table_vars, mode = "mixed", commas = TRUE, reliability_checks = FALSE, return_metrics = FALSE) { diff --git a/man/addint.Rd b/man/addint.Rd index 7ef6725..cc30b95 100644 --- a/man/addint.Rd +++ b/man/addint.Rd @@ -58,6 +58,15 @@ AP = RERI / RR11 S = (RR11 - 1) / ((RR10 - 1) + (RR01 - 1)) Confidence intervals for S are calculated on the log scale and then exponentiated. + +\strong{Scale caveat:} these measures are defined on the risk-ratio scale. +When the model supplies odds ratios (logistic) or hazard ratios (Cox), RERI, +AP, and S only approximate their risk-ratio counterparts if the outcome is +rare; for a common outcome the odds ratio overstates the risk ratio and the +measures can be biased. To target risk ratios directly, fit a log-binomial or +Poisson working model (see VanderWeele and Knol (2014) +). The delta-method Wald intervals here are +symmetric and approximate. } \examples{ \donttest{ diff --git a/man/addintlist.Rd b/man/addintlist.Rd index 5e33ef1..efbfb04 100644 --- a/man/addintlist.Rd +++ b/man/addintlist.Rd @@ -48,6 +48,9 @@ A `tibble` (data frame) summarizing the additive interaction results. \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. + \item `Scale`: The ratio scale the measures were computed on, inferred + from the model: "OR (logistic)", "HR (Cox)", "RR (log-binomial)", + "RR (Poisson)", or "ratio". } Returns an empty tibble if errors occur. } @@ -64,68 +67,39 @@ 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. + +RERI, AP, and the synergy index S are additive-interaction measures defined on +the \emph{risk-ratio} scale. When the model reports odds ratios (logistic) or +hazard ratios (Cox), these measures only approximate their risk-ratio +counterparts when the outcome is rare; for a common outcome the odds ratio +overstates the risk ratio and RERI/AP/S can be biased. The `Scale` column makes +the operating scale explicit. To estimate risk-ratio-based interaction +directly for a common outcome, fit a log-binomial or Poisson working model. } \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 - ) +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - # --- 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" - ) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # Print the results table - print(all_interactions_table, n=50) +# Saturated interaction model of sex and insulin use. +interaction_model <- survey::svyglm( + htn01 ~ sex * insulin + age, + design = design, family = quasibinomial() +) -} else { - print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.") -} +# RERI, AP and S for every factor-level combination. The Scale column flags +# that these are OR-based (see Details for the rare-outcome caveat). +addintlist( + model = interaction_model, + factor1_name = "sex", + factor2_name = "insulin", + measures = "all" +) } } \seealso{ diff --git a/man/jointeffects.Rd b/man/jointeffects.Rd index a5f702d..6d43b28 100644 --- a/man/jointeffects.Rd +++ b/man/jointeffects.Rd @@ -58,69 +58,25 @@ coefficients. } \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 ratio SE +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. Fit the Interaction Model --- - interaction_model_logit <- svyglm( - Hypertension_130 ~ Race1 * ObeseStatus + Age, - design = adult_design_binary, family = quasibinomial() - ) - - # --- 3. Run the function --- +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # Example 1: Output on Ratio scale (default) - joint_effects_ratio <- jointeffects( - interaction_model = interaction_model_logit, - factor1_name = "Race1", - factor2_name = "ObeseStatus", - scale = "ratio" - ) - print("--- Joint Effects (Ratio Scale) ---") - print(joint_effects_ratio, n = 50) +interaction_model <- survey::svyglm( + htn01 ~ sex * insulin + age, + design = design, family = quasibinomial() +) - # Example 2: Output on Log scale - joint_effects_log <- jointeffects( - interaction_model = interaction_model_logit, - factor1_name = "Race1", - factor2_name = "ObeseStatus", - scale = "log" - ) - print("--- Joint Effects (Log Scale) ---") - print(joint_effects_log, n = 50) -} +# Joint effects of sex x insulin on the ratio (OR) scale. +jointeffects( + interaction_model = interaction_model, + factor1_name = "sex", + factor2_name = "insulin", + scale = "ratio" +) } } diff --git a/man/plotint.Rd b/man/plotint.Rd index 86ddb83..739606b 100644 --- a/man/plotint.Rd +++ b/man/plotint.Rd @@ -101,94 +101,7 @@ plot_data <- plotint( legend_title = "log(Total Calories)" ) -# You can now access the plot data if you want it -# head(plot_data$plot_data) - - -# ================================================================= -# Example 2: Continuous x Continuous (NHANESraw) -# ================================================================= -library(NHANES) -library(dplyr) -library(survey) -library(emmeans) -library(ggplot2) - -# --- 3. DATA PREPARATION (NHANESraw) --- -data("NHANESraw", package = "NHANES") - -# Select relevant variables and filter for adults with complete cases -nhanes_clean <- NHANESraw \%>\% - filter(Age >= 20) \%>\% - mutate( - # Ensure factors are set correctly - Gender = factor(Gender), - Race1 = factor(Race1) - ) \%>\% - # Keep only the variables we need for this model - select( - BPSysAve, Age, BMI, Gender, Race1, - SDMVPSU, SDMVSTRA, WTMEC2YR - ) \%>\% - na.omit() # Listwise deletion for simplicity - -# --- 4. DESIGN & MODEL FITTING --- -# Create the survey design object -design_nhanes <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_clean -) - -# Fit the model: Interaction of Age * BMI on Geometric Mean of BP -# We use gaussian(link="log") to model Geometric Means -model_bp_interaction <- svyglm( - BPSysAve ~ Age * BMI + Gender + Race1, - design = design_nhanes, - family = gaussian(link = "log") -) - -# --- 5. PLOT THE INTERACTION --- - -# --- Plot 1: Standard Plot (using quantiles and CIs) --- -# This shows the effect of Age (x-axis) at 10th/50th/90th percentiles of BMI -print( - plotint( - model = model_bp_interaction, - effect = "Age", - moderator = "BMI", - data = nhanes_clean, - mod_scale = "quantile", # Use 10th/50th/90th percentiles for BMI - type = "response", # Plot the Geometric Mean (back-transform from log) - show_ci = TRUE, - bw = FALSE, - pval_position = "bottom.right", - xlab = "Age (Years)", - ylab = "Predicted Geometric Mean SBP (mmHg)", - legend_title = "BMI (Percentiles)" - ) -) - -# --- Plot 2: Cleaner Plot (using SD and CIs) --- -# This shows the effect of Age at Mean +/- 1 SD of BMI -# We turn off CIs for a cleaner look if they are too wide -print( - plotint( - model = model_bp_interaction, - effect = "Age", - moderator = "BMI", - data = nhanes_clean, - mod_scale = "sd", # Use Mean +/- 1 SD for BMI - type = "response", - show_ci = FALSE, # Turn off CI ribbons - bw = TRUE, # Use B&W linetypes - pval_position = "bottom.right", - xlab = "Age (Years)", - ylab = "Predicted Geometric Mean SBP (mmHg)", - legend_title = "BMI (Standard Deviations)" - ) -) +# You can now access the underlying plot data: +head(plot_data$plot_data) } } diff --git a/man/svydiag.Rd b/man/svydiag.Rd index 7bc1734..40f3e05 100644 --- a/man/svydiag.Rd +++ b/man/svydiag.Rd @@ -50,52 +50,17 @@ even if precisely estimated. It is better to rely on the standard error, p-value, and confidence interval width for reliability assessment. } \examples{ -# Ensure required packages are loaded -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE)) { +data(nhanes_mortality, package = "svyTable1") +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") - # 1. Prepare Data using the NHANES example - data(NHANESraw, package = "NHANES") - nhanes_adults_with_na <- NHANESraw \%>\% - dplyr::filter(Age >= 20) \%>\% - dplyr::mutate( - ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese")), - Race1 = factor(Race1) - ) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # Create a complete-case design object for the regression model - nhanes_complete <- nhanes_adults_with_na[complete.cases( - nhanes_adults_with_na[, c("ObeseStatus", "Age", "Race1")] - ), ] +fit <- survey::svyglm(htn01 ~ age + sex + race, + design = design, family = quasibinomial()) - adult_design_complete <- survey::svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_complete - ) - - # 2. Fit a survey-weighted logistic regression model - fit <- survey::svyglm( - ObeseStatus ~ Age + Race1, - design = adult_design_complete, - family = quasibinomial() - ) - - # 3. Get the reliability diagnostics table using the new function - diagnostics_table <- svydiag(fit) - - # Print the resulting table - print(diagnostics_table) - - # For a publication-ready table, pipe the result to kable() - if (requireNamespace("knitr", quietly = TRUE)) { - knitr::kable(diagnostics_table, - caption = "Reliability Diagnostics for NHANES Obesity Model", - digits = 3) - } -} +# Per-coefficient reliability diagnostics (SE, p-value, CI width, RSE). +svydiag(fit) } diff --git a/man/svytable1.Rd b/man/svytable1.Rd index f0754cb..f810feb 100644 --- a/man/svytable1.Rd +++ b/man/svytable1.Rd @@ -71,50 +71,29 @@ optional comma formatting, and optional reliability checks based on NCHS Data Presentation Standards. } \examples{ -# Ensure required packages are loaded -if (requireNamespace("survey", quietly = TRUE) && - requireNamespace("NHANES", quietly = TRUE) && - requireNamespace("dplyr", quietly = TRUE)) { +data(nhanes_mortality, package = "svyTable1") - # 1. Prepare Data using NHANES - data(NHANESraw, package = "NHANES") - nhanes_adults_with_na <- NHANESraw \%>\% - dplyr::filter(Age >= 20) \%>\% - dplyr::mutate( - ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese")) - ) - - adult_design_with_na <- survey::svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_adults_with_na - ) - - # 2. Basic Example: Create a simple Table 1 - vars_to_summarize <- c("Age", "Gender", "Race1", "Education") - table1 <- svytable1( - design = adult_design_with_na, - strata_var = "ObeseStatus", - table_vars = vars_to_summarize - ) - print(table1) - - # 3. Advanced Example: Use reliability checks and get detailed metrics - results <- svytable1( - design = adult_design_with_na, - strata_var = "ObeseStatus", - table_vars = vars_to_summarize, - reliability_checks = TRUE, - return_metrics = TRUE - ) +design <- survey::svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality +) - # View the table with unreliable estimates suppressed (*) - print(results$formatted_table) +# Basic survey-weighted Table 1, stratified by hypertension status. +svytable1( + design = design, + strata_var = "htn", + table_vars = c("age", "sex", "race", "smoking") +) - # View the detailed report card of reliability checks - print(results$reliability_metrics) -} +# With NCHS reliability checks plus a detailed metrics table. Estimates that +# fail the reliability rules are shown as "*". +results <- svytable1( + design = design, + strata_var = "htn", + table_vars = c("age", "sex", "race", "smoking"), + reliability_checks = TRUE, + return_metrics = TRUE +) +results$formatted_table +head(results$reliability_metrics) } diff --git a/tests/testthat/test-interaction.R b/tests/testthat/test-interaction.R index 7a9fad4..6501613 100644 --- a/tests/testthat/test-interaction.R +++ b/tests/testthat/test-interaction.R @@ -17,6 +17,13 @@ test_that("addintlist returns the three additive-interaction measures", { expect_setequal(unique(al$Measure), c("RERI", "AP", "S")) }) +test_that("addintlist labels the effect-measure scale (OR for logistic)", { + al <- addintlist(model = st1_int_fit(), factor1_name = "sex", + factor2_name = "insulin", measures = "all") + expect_true("Scale" %in% names(al)) + expect_equal(unique(al$Scale), "OR (logistic)") +}) + test_that("reportint(output = 'list') returns the documented panels", { fit <- st1_int_fit() suppressMessages(invisible(utils::capture.output( diff --git a/vignettes/v1-svyTable1-basics.Rmd b/vignettes/v1-svyTable1-basics.Rmd index 61a141d..67f54e5 100644 --- a/vignettes/v1-svyTable1-basics.Rmd +++ b/vignettes/v1-svyTable1-basics.Rmd @@ -9,7 +9,12 @@ vignette: > %\VignetteEncoding{UTF-8} --- -`{r setup} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )` +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` The `svyTable1` package provides a suite of functions for creating publication-ready tables from complex analytical results. This vignette @@ -33,42 +38,39 @@ shown for categorical variables: This structure prevents common reporting problems such as showing weighted counts (which do not represent real sample sizes) or hiding sparse categories that contribute instability. -For continuous variables, the function reports weighted means, standard errors, and confidence intervals. It also offers optional reliability checks derived from National Center for Health Statistics (NCHS) guidance, including standard error and relative standard error thresholds. +For continuous variables, the function reports weighted means and standard deviations. It also offers optional reliability checks derived from National Center for Health Statistics (NCHS) guidance. -## Example: Using NHANES Data (2009--2012) +## Example: The Bundled `nhanes_mortality` Cohort -The examples below use data from the NHANES raw dataset to demonstrate how to prepare a survey design and generate descriptive tables. The code illustrates both incomplete data situations and complete case summaries. +The examples below use the bundled `nhanes_mortality` dataset (a processed NHANES cohort of diabetic adults). The code illustrates both incomplete-data situations and complete-case summaries, stratified by hypertension status (`htn`). ### 1. Data Preparation ```{r packages, message = FALSE, warning = FALSE} library(svyTable1) library(survey) -library(dplyr) -library(NHANES) library(knitr) -library(tidyr) + +options(survey.lonely.psu = "adjust") ``` -``` {r data} -data(NHANESraw) +```{r data} +data(nhanes_mortality) -nhanes_adults_with_na <- NHANESraw %>% - dplyr::filter(Age >= 20) %>% - mutate( - ObeseStatus = factor( - ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese") - ) - ) +# This dataset is complete; induce some missingness to illustrate how +# svytable1() reports it. +set.seed(42) +nhanes_miss <- nhanes_mortality +nhanes_miss$smoking[sample(nrow(nhanes_miss), 150)] <- NA +nhanes_miss$bmi.cat[sample(nrow(nhanes_miss), 100)] <- NA -adult_design_with_na <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, +design_miss <- svydesign( + id = ~psu, + strata = ~strata, + weights = ~survey_weight, nest = TRUE, - data = nhanes_adults_with_na + data = nhanes_miss ) ``` @@ -78,15 +80,12 @@ adult_design_with_na <- svydesign( Variables in large health surveys often contain missing values. `svytable1()` automatically identifies missingness and calculates unweighted Ns that reflect real denominators. Weighted percentages are based only on cases with observed data, consistent with standard survey practice. -``` {r miss} -vars_with_missing <- c( - "Age", "Gender", "Race1", "Education", "HHIncome", - "TotChol", "SleepHrsNight", "SmokeNow" -) +```{r miss} +vars_with_missing <- c("age", "sex", "race", "smoking", "bmi.cat") table_with_missing <- svytable1( - design = adult_design_with_na, - strata_var = "ObeseStatus", + design = design_miss, + strata_var = "htn", table_vars = vars_with_missing ) @@ -98,28 +97,21 @@ knitr::kable( #### Example B: Summarizing Complete Data -If you prefer to restrict the analysis to respondents with complete data for a chosen set of variables, you can construct a complete analytic dataset and supply it to a new `svydesign` object. The function then calculates summary statistics without any missing categories. +If you prefer to restrict the analysis to respondents with complete data, you can supply a design built on the complete dataset. The function then calculates summary statistics without any missing categories. -``` {r complete} -vars_for_complete_table <- c( - "Age", "Gender", "Race1", "BPSysAve", "Pulse", "BMI" -) - -nhanes_adults_complete <- nhanes_adults_with_na %>% - drop_na(all_of(vars_for_complete_table)) - -adult_design_complete <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, +```{r complete} +design_complete <- svydesign( + id = ~psu, + strata = ~strata, + weights = ~survey_weight, nest = TRUE, - data = nhanes_adults_complete + data = nhanes_mortality ) table_without_missing <- svytable1( - design = adult_design_complete, - strata_var = "ObeseStatus", - table_vars = c("Age", "Gender", "Race1", "BPSysAve", "Pulse") + design = design_complete, + strata_var = "htn", + table_vars = c("age", "sex", "race") ) knitr::kable( @@ -128,20 +120,20 @@ knitr::kable( ) ``` -### 3. Checking Estimate Reliability for Proportions (NCHS Standards) +### 3. Checking Estimate Reliability (NCHS Standards) Federal statistical agencies discourage reporting proportions that are unstable or based on very small samples. When you request `reliability_checks = TRUE`, `svytable1()` returns: -- The formatted table with suppressed estimates. +- The formatted table with unreliable estimates suppressed (shown as `*`). - A companion metrics table showing which estimates failed reliability standards. -These checks follow NCHS guidance for the evaluation of proportions in complex surveys. Estimates with relative standard error values above recommended thresholds are marked or suppressed. +These checks follow NCHS guidance for the evaluation of proportions in complex surveys [@nchs_reliability_estimates]. -``` {r rel} +```{r rel} results_list <- svytable1( - design = adult_design_with_na, - strata_var = "ObeseStatus", + design = design_miss, + strata_var = "htn", table_vars = vars_with_missing, reliability_checks = TRUE, return_metrics = TRUE @@ -151,26 +143,27 @@ knitr::kable( results_list$formatted_table, caption = "Table 3: Table with NCHS Reliability Checks Applied" ) - -knitr::kable( - results_list$reliability_metrics[ - results_list$reliability_metrics$suppressed == TRUE, ], - caption = "Reliability Metrics for Suppressed Estimates" -) ``` -### 4. Checking Estimate Reliability for Means (NCHS Standards) - -For numeric variables, `svytable1()` performs an additional check: - -- `fail_rse_30`: Flags means whose relative standard error (RSE) is at - least 30 percent. +The companion metrics table records, for every evaluated cell, the unweighted +*n*, design degrees of freedom, design effect, effective sample size, +confidence-interval bounds, and the specific reliability rules that were +triggered. Filtering it to the suppressed cells explains each `*`: -RSE is calculated as: +```{r relmetrics} +suppressed <- results_list$reliability_metrics[ + results_list$reliability_metrics$suppressed == TRUE, ] - RSE = (Standard Error / Estimate) * 100 +knitr::kable( + head(suppressed, 10), + caption = "Reliability Metrics for Suppressed Estimates", + digits = 2, + row.names = FALSE +) +``` -An RSE ≥ 30 percent indicates an unreliable estimate based on NCHS -guidance. +For numeric variables, the relevant rule is the relative standard error +(RSE = standard error / estimate); an estimate is flagged when its RSE reaches +30 percent, following NCHS guidance. ## References diff --git a/vignettes/v2-model-diagnostics.Rmd b/vignettes/v2-model-diagnostics.Rmd index 76d5ee1..91ae08b 100644 --- a/vignettes/v2-model-diagnostics.Rmd +++ b/vignettes/v2-model-diagnostics.Rmd @@ -9,7 +9,12 @@ vignette: > %\VignetteEncoding{UTF-8} --- -`{r setup} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )` +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` This vignette covers functions in the `svyTable1` package for diagnosing fitted survey regression models: @@ -23,50 +28,34 @@ These tools provide a structured way to assess whether model coefficients are pr ## 1. Data and Model Preparation -The examples use complete-case NHANES data from adults aged 20 years or older. A survey-weighted logistic regression model is fitted to predict obesity. This model is used throughout the vignette to illustrate each diagnostic function. +The examples use the bundled `nhanes_mortality` dataset (a processed NHANES +cohort of diabetic adults). A survey-weighted logistic regression model is +fitted to predict hypertension (`htn`). This model is used throughout the +vignette to illustrate each diagnostic function. ```{r packages, warning = FALSE, message = FALSE} library(svyTable1) library(survey) -library(dplyr) -library(NHANES) library(knitr) -library(tidyr) library(ggplot2) ``` -``` {r data} -data(NHANESraw) - -nhanes_adults_with_na <- NHANESraw %>% - filter(Age >= 20) %>% - mutate( - ObeseStatus = factor( - ifelse(BMI >= 30, "Obese", "Not Obese"), - levels = c("Not Obese", "Obese") - ) - ) - -vars_for_complete_table <- c( - "Age", "Gender", "Race1", "BPSysAve", - "Pulse", "BMI", "ObeseStatus" -) - -nhanes_adults_complete <- nhanes_adults_with_na %>% - drop_na(all_of(vars_for_complete_table)) +```{r data} +data(nhanes_mortality) +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") -adult_design_complete <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, +adult_design <- svydesign( + id = ~psu, + strata = ~strata, + weights = ~survey_weight, nest = TRUE, - data = nhanes_adults_complete + data = nhanes_mortality ) -fit_obesity <- svyglm( - ObeseStatus ~ Age + Gender + Race1, - design = adult_design_complete, +fit_htn <- svyglm( + htn01 ~ age + sex + smoking, + design = adult_design, family = quasibinomial() ) ``` @@ -76,12 +65,12 @@ fit_obesity <- svyglm( `svydiag()` produces diagnostics for regression coefficients, including p-values, precision measures, and relative standard errors. -``` {r diag} -diagnostics_table <- svydiag(fit_obesity) +```{r diag} +diagnostics_table <- svydiag(fit_htn) knitr::kable( diagnostics_table, - caption = "Reliability Diagnostics for NHANES Obesity Model", + caption = "Reliability Diagnostics for the Hypertension Model", digits = 3 ) ``` @@ -98,48 +87,48 @@ knitr::kable( `svygof()` implements the Archer and Lemeshow goodness of fit procedure for complex survey logistic regression [@archer2006goodness]. This test adapts the Hosmer and Lemeshow approach to designs with weights, strata, and clustering. -The function groups respondents by predicted risk and compares observed and expected outcomes while incorporating survey weights. A small p-value provides evidence that predicted risks diverge from observed risks in one or more groups. +The function groups respondents by predicted risk and compares observed and expected outcomes while incorporating survey weights. A small p-value provides evidence that predicted risks diverge from observed risks in one or more groups. Note that `svygof()` expects a standard `svydesign` object (not a replicate-weights design). This test helps identify lack of fit that may result from omitted interactions or nonlinearities in the model. -``` {r gof} -gof_results <- svygof(fit_obesity, adult_design_complete) +```{r gof} +gof_results <- svygof(fit_htn, adult_design) knitr::kable( gof_results, - caption = "Archer-Lemeshow Goodness-of-Fit Test for Obesity Model" + caption = "Archer-Lemeshow Goodness-of-Fit Test for the Hypertension Model" ) ``` -A significant p-value indicates poor model fit across predicted risk -groups. +A non-significant p-value is consistent with adequate fit across predicted +risk groups; it does not prove the model is correctly specified. ## 4. Model Performance (`svyAUC()`) -`svyAUC()` estimates the area under the ROC curve using replicate-weight designs. This produces an AUC estimate and confidence interval that appropriately reflect survey sampling. +`svyAUC()` estimates the area under the ROC curve using replicate-weight designs. This produces an AUC estimate and confidence interval that appropriately reflect survey sampling. The confidence interval is built on the logit scale, so it always stays within `[0, 1]`. -The function also returns data for constructing ROC curves. A plot of the ROC curve, along with the AUC estimate, demonstrates how well the model discriminates between obese and non-obese respondents. +The function also returns data for constructing ROC curves. A plot of the ROC curve, along with the AUC estimate, demonstrates how well the model discriminates between respondents with and without hypertension. A value close to 0.5 indicates little discriminatory value, while higher values indicate stronger predictive performance. -``` {r auc} -rep_design <- as.svrepdesign(adult_design_complete) +```{r auc} +rep_design <- as.svrepdesign(adult_design) -fit_obesity_rep <- svyglm( - ObeseStatus ~ Age + Gender + Race1, +fit_htn_rep <- svyglm( + htn01 ~ age + sex + smoking, design = rep_design, family = quasibinomial() ) -auc_results_list <- svyAUC(fit_obesity_rep, rep_design, plot = TRUE) +auc_results_list <- svyAUC(fit_htn_rep, rep_design, plot = TRUE) knitr::kable( auc_results_list$summary, - caption = "Design-Correct AUC for Obesity Model" + caption = "Design-Correct AUC for the Hypertension Model" ) ggplot(auc_results_list$roc_data, aes(x = FPR, y = TPR)) + - geom_line(linewidth = 1) + # Changed 'size' to 'linewidth' + geom_line(linewidth = 1) + geom_abline(linetype = "dashed") + labs( title = "Survey-Weighted ROC Curve", diff --git a/vignettes/v3-interaction-analysis.Rmd b/vignettes/v3-interaction-analysis.Rmd index 8f32578..dea3d4e 100644 --- a/vignettes/v3-interaction-analysis.Rmd +++ b/vignettes/v3-interaction-analysis.Rmd @@ -23,13 +23,13 @@ This vignette demonstrates how to analyze and present interaction effects using - Generating automated multi-panel reports (`reportint()`) - Visualizing interactions involving continuous variables (`plotint()`) -We will use NHANES data to explore these functions. +We use the bundled `nhanes_mortality` dataset (a processed NHANES cohort of diabetic adults) throughout. ## Analyzing Interaction Effects (Categorical) The `svyTable1` package includes several functions to help analyze and present interaction effects between two categorical variables in survey-weighted regression models (`svyglm` or `svycoxph`). These functions facilitate demonstrating the equivalence of different modeling approaches and calculating measures of additive interaction. -### Example: Interaction between Race and Obesity on Hypertension (NHANES) +### Example: Interaction between Sex and Insulin Use on Hypertension For two categorical variables, `svyTable1` provides functions that summarize joint, simple, and additive interaction effects from survey-weighted logistic or Cox models. These summaries help demonstrate the connection between the saturated interaction model and the equivalent joint-indicator model. @@ -38,64 +38,50 @@ For two categorical variables, `svyTable1` provides functions that summarize joi ```{r packages, warning = FALSE, message = FALSE} library(svyTable1) library(survey) -library(dplyr) -library(NHANES) library(knitr) -library(tidyr) library(emmeans) library(ggplot2) + +options(survey.lonely.psu = "adjust") ``` -The example uses NHANES data to examine whether the association between obesity and hypertension differs by race. A survey-weighted logistic model with an interaction term is fitted, and a matched joint-indicator model is created to demonstrate equivalence. +The example examines whether the association between insulin use and hypertension differs by sex. A survey-weighted logistic model with an interaction term is fitted, and a matched joint-indicator model is created to demonstrate equivalence. ```{r data} -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, +data(nhanes_mortality) +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") + +design <- svydesign( + id = ~psu, + strata = ~strata, + weights = ~survey_weight, nest = TRUE, - data = nhanes_adults_processed + data = nhanes_mortality ) +# Saturated interaction model interaction_model_logit <- svyglm( - Hypertension_130 ~ Race1 * ObeseStatus + Age, - design = adult_design_binary, + htn01 ~ sex * insulin + age, + design = design, family = quasibinomial() ) -adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = interaction(Race1, ObeseStatus, sep = "_", drop = TRUE) +# Equivalent joint-indicator model +nhanes_mortality$sex_insulin <- interaction( + nhanes_mortality$sex, nhanes_mortality$insulin, sep = "_" ) - -adult_design_binary <- update(adult_design_binary, - Race1_ObeseStatus = relevel(Race1_ObeseStatus, ref = "White_Not Obese") +design_joint <- svydesign( + id = ~psu, strata = ~strata, weights = ~survey_weight, + nest = TRUE, data = nhanes_mortality ) - joint_model_logit <- svyglm( - Hypertension_130 ~ Race1_ObeseStatus + Age, - design = adult_design_binary, + htn01 ~ sex_insulin + age, + design = design_joint, family = quasibinomial() ) -f1_levels <- levels(adult_design_binary$variables$Race1) -f2_levels <- levels(adult_design_binary$variables$ObeseStatus) +f1_levels <- c("Male", "Female") # reference level first +f2_levels <- c("No", "Yes") # reference level first ``` ## 2. Additive Interaction Measures @@ -108,19 +94,19 @@ f2_levels <- levels(adult_design_binary$variables$ObeseStatus) - S (synergy index) -These measures quantify departure from additivity on the ratio scale and help assess whether the joint effect of two exposures exceeds what would be expected from independent effects. +These measures quantify departure from additivity. Note that for logistic models they are computed from odds ratios, which approximate risk ratios only when the outcome is rare; interpret them with care for common outcomes. ```{r add} all_interactions_table <- addintlist( model = interaction_model_logit, - factor1_name = "Race1", - factor2_name = "ObeseStatus", + factor1_name = "sex", + factor2_name = "insulin", measures = "all", digits = 3 ) knitr::kable(all_interactions_table, - caption = "Additive Interaction Measures (RERI, AP, S) for Race x Obesity on Hypertension", + caption = "Additive Interaction Measures (RERI, AP, S) for Sex x Insulin on Hypertension", digits = 3) ``` @@ -131,37 +117,26 @@ knitr::kable(all_interactions_table, ```{r joint} joint_effects_ratio <- jointeffects( interaction_model = interaction_model_logit, - factor1_name = "Race1", - factor2_name = "ObeseStatus", + factor1_name = "sex", + factor2_name = "insulin", scale = "ratio", digits = 2) knitr::kable(joint_effects_ratio, caption = "Joint Effects (OR Scale) Estimated from Interaction Model", digits = 2) - -joint_effects_log <- jointeffects( - interaction_model = interaction_model_logit, - factor1_name = "Race1", - factor2_name = "ObeseStatus", - scale = "log", - digits = 3) - -knitr::kable(joint_effects_log, - caption = "Joint Effects (Log-OR Scale) Estimated from Interaction Model", - digits = 3) ``` ## 4. Simple Effects -`inteffects()` provides stratum-specific effects. For example, it reports the effect of obesity within each racial group and the effect of race within each obesity category. This helps identify which group-specific comparisons drive the interaction. +`inteffects()` provides stratum-specific effects: the effect of insulin within each sex, and the effect of sex within each insulin category. This helps identify which group-specific comparisons drive the interaction. ```{r simpeff} simple_effects_ratio <- inteffects( joint_model = joint_model_logit, - joint_var_name = "Race1_ObeseStatus", - factor1_name = "Race1", - factor2_name = "ObeseStatus", + joint_var_name = "sex_insulin", + factor1_name = "sex", + factor2_name = "insulin", factor1_levels = f1_levels, factor2_levels = f2_levels, level_separator = "_", @@ -171,58 +146,26 @@ simple_effects_ratio <- inteffects( knitr::kable(simple_effects_ratio, caption = "Simple Effects (OR Scale) Estimated from Joint Model", digits = 2) - -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", - digits = 3) - -knitr::kable(simple_effects_log, - caption = "Simple Effects (Log-OR Scale) Estimated from Joint Model", - digits = 3) ``` ## 5. Automated Multi-Panel Report -`reportint()` compiles the results into a structured multi-panel summary that includes: - -- joint effects -- simple effects -- additive interaction measures -- multiplicative interaction estimates -- model details +`reportint()` compiles the results into a structured multi-panel summary that includes joint effects, simple effects, additive interaction measures, multiplicative interaction estimates, and model details. ```{r auto} interaction_report <- reportint( interaction_model = interaction_model_logit, + factor1_name = "sex", + factor2_name = "insulin", output = "list" ) knitr::kable(interaction_report$joint_effects, caption = "Panel A (Joint Effects) from reportint()") -knitr::kable(interaction_report$stratum_specific_effects, - caption = "Panel B (Simple Effects) from reportint()") - knitr::kable(interaction_report$additive_interaction, caption = "Panel C (Additive Interaction) from reportint()") -knitr::kable(interaction_report$model_details, - caption = "Panel D (Model Details) from reportint()") - -knitr::kable(interaction_report$multiplicative_scale, - caption = "Multiplicative Interaction (RORs) from reportint()") - -knitr::kable(interaction_report$effect_modification_report, - caption = "Publication Report: Effect Modification Analysis", - align = 'l') - knitr::kable(interaction_report$Interaction_report, caption = "Publication Report: Full Interaction Analysis", align = 'l') @@ -232,51 +175,28 @@ This panel format follows reporting recommendations [@knol2012recommendations]. ## Visualizing Interaction Effects with `plotint` -`plotint()` produces plots of predicted outcomes across values of a continuous predictor for selected levels of a continuous or categorical moderator. The function uses survey-weighted marginal estimates, supports both model and data driven moderator scales, and returns data and a plot object. +`plotint()` produces plots of predicted outcomes across values of a continuous predictor for selected levels of a continuous or categorical moderator. The function uses survey-weighted marginal estimates, supports both model- and data-driven moderator scales, and returns the underlying plot data and a plot object. -This approach is helpful when evaluating interactions between continuous variables, where reporting all numerical contrasts is often impractical. +This approach is helpful when evaluating interactions between continuous variables, where reporting all numerical contrasts is impractical. Here we model total caloric intake as a function of an age-by-carbohydrate interaction. ```{r viz} -data("NHANESraw") -nhanes_clean <- NHANESraw %>% - filter(Age >= 20) %>% - mutate( - Gender = factor(Gender), - Race1 = factor(Race1) - ) %>% - select( - BPSysAve, Age, BMI, Gender, Race1, - SDMVPSU, SDMVSTRA, WTMEC2YR - ) %>% - na.omit() - -design_nhanes <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, - nest = TRUE, - data = nhanes_clean -) - -model_bp_interaction <- svyglm( - BPSysAve ~ Age * BMI + Gender + Race1, - design = design_nhanes, - family = gaussian(link = "log") +model_cal_interaction <- svyglm( + cal.total ~ age * carbohyd + sex, + design = design, + family = gaussian() ) plot_data_list <- plotint( - model = model_bp_interaction, - effect = "Age", - moderator = "BMI", - data = nhanes_clean, + model = model_cal_interaction, + effect = "age", + moderator = "carbohyd", + data = nhanes_mortality, mod_scale = "quantile", type = "response", show_ci = TRUE, - bw = FALSE, - pval_position = "bottom.right", - xlab = "Age (Years)", - ylab = "Predicted Geometric Mean SBP (mmHg)", - legend_title = "BMI (Percentiles)" + xlab = "Age (years)", + ylab = "Predicted total caloric intake (kcal)", + legend_title = "Carbohydrate (percentiles)" ) head(plot_data_list$plot_data) diff --git a/vignettes/v4-survival-analysis.Rmd b/vignettes/v4-survival-analysis.Rmd index 4a63d44..434adc8 100644 --- a/vignettes/v4-survival-analysis.Rmd +++ b/vignettes/v4-survival-analysis.Rmd @@ -33,8 +33,9 @@ library(knitr) library(survival) library(ggplot2) library(dplyr) -library(mitools) library(mice) + +options(survey.lonely.psu = "adjust") ``` `svykmplot()` produces Kaplan–Meier survival estimates that correctly incorporate survey weights, strata, and clusters. The function returns both a plot and a number-at-risk table. The plot supports automatic or custom colors, p-value display, and censoring indicators. @@ -223,7 +224,10 @@ svycoxph_CE( - Multiple Imputation Case (`svycoxph_CE_mi`) -```{r mi} +The multiple-imputation pipeline relies on `mitools::imputationList()`, so this +chunk only runs when the suggested **mitools** package is installed. + +```{r mi, eval = requireNamespace("mitools", quietly = TRUE)} data(nhanes_mortality) analytic_design <- svydesign( @@ -276,7 +280,7 @@ imputed_data <- mice( impdata_long <- mice::complete(imputed_data, "long", include = FALSE) -allImputations_main <- imputationList(split(impdata_long, impdata_long$.imp)) +allImputations_main <- mitools::imputationList(split(impdata_long, impdata_long$.imp)) design_main <- svydesign( strata = ~strata, id = ~psu, @@ -300,7 +304,7 @@ impdata_long_split <- survSplit(Surv(stime, status) ~ ., episode = "tgroup", id = "split_id") -allImputations_split <- imputationList(split(impdata_long_split, impdata_long_split$.imp)) +allImputations_split <- mitools::imputationList(split(impdata_long_split, impdata_long_split$.imp)) design_split <- svydesign( strata = ~strata, diff --git a/vignettes/v5-multiple-imputation.Rmd b/vignettes/v5-multiple-imputation.Rmd index f119e72..ff96e93 100644 --- a/vignettes/v5-multiple-imputation.Rmd +++ b/vignettes/v5-multiple-imputation.Rmd @@ -22,57 +22,60 @@ This vignette demonstrates how to use the **svypooled()** function to format poo Beyond survey data, a common analytical challenge is presenting results from models fitted to multiply imputed datasets. The **svyTable1** package provides the `svypooled()` function to address this. It takes a pooled model object (`mipo`) from **mice** and produces publication-ready regression tables. -A key feature of `svypooled()` is its ability to create a fallacy-safe table. This approach displays results for the primary exposure only while listing adjustment covariates in a footnote, preventing misinterpretation of covariate p-values. +A key feature of `svypooled()` is its ability to create a fallacy-safe table. This approach displays results for the primary exposure only while listing adjustment covariates in a footnote, preventing misinterpretation of covariate p-values [@westreich2013table]. ## Example: Handling Missing Data in a Survey Design -The following sections illustrate the recommended workflow. Data are taken from the NHANES example dataset supplied in the `NHANES` package. +The following sections illustrate the recommended workflow using the bundled +`nhanes_mortality` dataset (a processed NHANES cohort of diabetic adults). ### Data Preparation ```{r packages, warning = FALSE, message = FALSE} library(svyTable1) library(mice) -library(dplyr) library(survey) -library(NHANES) library(knitr) -library(mitools) + +options(survey.lonely.psu = "adjust") ``` -We prepare an analytic dataset, define obesity, categorize age, define a binary smoking variable, and select variables needed for both imputation and the survey design. +We build an analytic dataset, define a binary hypertension outcome, and -- so +the imputation step has something to do -- induce some missingness in a +covariate. ```{r data} -data(NHANESraw, package = "NHANES") -nhanes_analytic <- NHANESraw %>% - filter(Age >= 20 & WTMEC2YR > 0) %>% - mutate( - Obese = factor(ifelse(BMI >= 30, "Yes", "No"), levels = c("No", "Yes")), - 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, AgeCat, Smoke100, Education, SDMVPSU, SDMVSTRA, WTMEC2YR) +data(nhanes_mortality) +nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes") + +analytic <- nhanes_mortality[, c("htn01", "sex", "age", "smoking", + "psu", "strata", "survey_weight")] + +set.seed(123) +analytic$age[sample(nrow(analytic), 200)] <- NA ``` ## Pooling from Imputed Data -We run `mice()`, fit a survey-weighted logistic regression model to each imputed dataset using `svyglm()`, and combine the results using `pool()`. This pooled model object is accepted directly by `svypooled()`. +We run `mice()`, fit a survey-weighted logistic regression model to each imputed dataset using `svyglm()`, and combine the results using `pool()` (Rubin's rules). This pooled model object is accepted directly by `svypooled()`. + +A small number of imputations is used here only to keep the vignette fast; for a real analysis use `m >= 20` and check convergence. ```{r pool} -imputed_data <- mice(nhanes_analytic, m = 2, maxit = 2, seed = 123, printFlag = FALSE) +imputed_data <- mice(analytic, m = 2, maxit = 2, seed = 123, printFlag = FALSE) fit_list <- list() -for (i in 1:imputed_data$m) { +for (i in seq_len(imputed_data$m)) { completed_data <- complete(imputed_data, i) design_i <- svydesign( - id = ~SDMVPSU, - strata = ~SDMVSTRA, - weights = ~WTMEC2YR, + id = ~psu, + strata = ~strata, + weights = ~survey_weight, nest = TRUE, data = completed_data ) fit_list[[i]] <- svyglm( - Obese ~ Smoke100 + AgeCat + Education, + htn01 ~ sex + age + smoking, design = design_i, family = quasibinomial() ) @@ -83,17 +86,17 @@ pooled_results <- pool(fit_list) ## Generating Tables with `svypooled()` -### Fallacy-Safe Table +### Fallacy-Safe Table This is the recommended format [@westreich2013table]. It presents results only for the main exposure and lists adjustment variables as a single footnote. ```{r t2} svypooled( pooled_model = pooled_results, - main_exposure = "Smoke100", - adj_var_names = c("AgeCat", "Education"), + main_exposure = "sex", + adj_var_names = c("age", "smoking"), measure = "OR", - title = "Adjusted Odds of Obesity (Fallacy-Safe)" + title = "Adjusted Odds of Hypertension (Fallacy-Safe)" ) ``` @@ -104,10 +107,10 @@ This version displays results for every variable in the model. It is useful for ```{r t2pool} svypooled( pooled_model = pooled_results, - main_exposure = "Smoke100", - adj_var_names = c("AgeCat", "Education"), + main_exposure = "sex", + adj_var_names = c("age", "smoking"), measure = "OR", - title = "Adjusted Odds of Obesity (Full Table for Appendix)", + title = "Adjusted Odds of Hypertension (Full Table for Appendix)", fallacy_safe = FALSE ) ```