Self-contain vignettes/examples on bundled data; label interaction scale

DESCRIPTION

@@ -57,7 +57,6 @@ Imports:
     WeightedROC,
     broom
 Suggests:
-    NHANES,
     mitools,
     testthat (>= 3.0.0)
 VignetteBuilder: knitr

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)
+#' <doi:10.1515/em-2013-0005>). 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.

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"
+}

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,

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,