Skip to content

Commit 0c5f4f4

Browse files
ehsanxehsanx
authored
Merge pull request #2 from ehsanx/improve/vignettes-and-scale
Self-contain vignettes/examples on bundled data; label interaction scale
2 parents 0aab7e3 + 5ff07e5 commit 0c5f4f4

19 files changed

Lines changed: 401 additions & 874 deletions

‎DESCRIPTION‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,6 @@ Imports:
5757
WeightedROC,
5858
broom
5959
Suggests:
60-
NHANES,
6160
mitools,
6261
testthat (>= 3.0.0)
6362
VignetteBuilder: knitr

‎R/addint.R‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,15 @@
2222
#'
2323
#' Confidence intervals for S are calculated on the log scale and then exponentiated.
2424
#'
25+
#' \strong{Scale caveat:} these measures are defined on the risk-ratio scale.
26+
#' When the model supplies odds ratios (logistic) or hazard ratios (Cox), RERI,
27+
#' AP, and S only approximate their risk-ratio counterparts if the outcome is
28+
#' rare; for a common outcome the odds ratio overstates the risk ratio and the
29+
#' measures can be biased. To target risk ratios directly, fit a log-binomial or
30+
#' Poisson working model (see VanderWeele and Knol (2014)
31+
#' <doi:10.1515/em-2013-0005>). The delta-method Wald intervals here are
32+
#' symmetric and approximate.
33+
#'
2534
#' @param model A fitted model object (e.g., `svycoxph`, `svyglm`).
2635
#' @param type Character string: `"joint"` if using a combined categorical variable,
2736
#' `"interaction"` if using main effects and a product term.

‎R/addintlist.R‎

Lines changed: 52 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -40,9 +40,21 @@
4040
#' \item `SE`: Standard Error (Note: For S, this is SE of log(S)).
4141
#' \item `CI_low`: Lower confidence interval bound.
4242
#' \item `CI_upp`: Upper confidence interval bound.
43+
#' \item `Scale`: The ratio scale the measures were computed on, inferred
44+
#' from the model: "OR (logistic)", "HR (Cox)", "RR (log-binomial)",
45+
#' "RR (Poisson)", or "ratio".
4346
#' }
4447
#' Returns an empty tibble if errors occur.
4548
#'
49+
#' @details
50+
#' RERI, AP, and the synergy index S are additive-interaction measures defined on
51+
#' the \emph{risk-ratio} scale. When the model reports odds ratios (logistic) or
52+
#' hazard ratios (Cox), these measures only approximate their risk-ratio
53+
#' counterparts when the outcome is rare; for a common outcome the odds ratio
54+
#' overstates the risk ratio and RERI/AP/S can be biased. The `Scale` column makes
55+
#' the operating scale explicit. To estimate risk-ratio-based interaction
56+
#' directly for a common outcome, fit a log-binomial or Poisson working model.
57+
#'
4658
#' @seealso \code{\link{addint}}
4759
#'
4860
#' @importFrom stats coef
@@ -55,65 +67,28 @@
5567
#'
5668
#' @examples
5769
#' \donttest{
58-
#' # --- Load required libraries for the example ---
59-
#' # Ensure the 'addint' function is defined or loaded from the package
60-
#' # source("R/addint.R") # If running interactively before package build
61-
#'
62-
#' if (requireNamespace("survey", quietly = TRUE) &&
63-
#' requireNamespace("NHANES", quietly = TRUE) &&
64-
#' requireNamespace("dplyr", quietly = TRUE) &&
65-
#' requireNamespace("tidyr", quietly = TRUE) &&
66-
#' requireNamespace("msm", quietly = TRUE)) {
67-
#'
68-
#' library(survey)
69-
#' library(NHANES)
70-
#' library(dplyr)
71-
#' library(tidyr)
72-
#' library(msm)
73-
#'
74-
#' # --- 1. Data Preparation (NHANES Example) ---
75-
#' data(NHANESraw)
76-
#'
77-
#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130",
78-
#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR")
79-
#'
80-
#' nhanes_adults_processed <- NHANESraw %>%
81-
#' filter(Age >= 20) %>%
82-
#' mutate(
83-
#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"),
84-
#' levels = c("Not Obese", "Obese")),
85-
#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"),
86-
#' levels = c("No", "Yes")),
87-
#' Race1 = relevel(as.factor(Race1), ref = "White")
88-
#' ) %>%
89-
#' select(all_of(vars_needed)) %>%
90-
#' drop_na()
70+
#' data(nhanes_mortality, package = "svyTable1")
71+
#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes")
9172
#'
92-
#' adult_design_binary <- svydesign(
93-
#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR,
94-
#' nest = TRUE, data = nhanes_adults_processed
95-
#' )
73+
#' design <- survey::svydesign(
74+
#' id = ~psu, strata = ~strata, weights = ~survey_weight,
75+
#' nest = TRUE, data = nhanes_mortality
76+
#' )
9677
#'
97-
#' # --- 2. Fit Interaction Term Model ---
98-
#' interaction_model_logit <- svyglm(
99-
#' Hypertension_130 ~ Race1 * ObeseStatus + Age,
100-
#' design = adult_design_binary, family = quasibinomial()
101-
#' )
78+
#' # Saturated interaction model of sex and insulin use.
79+
#' interaction_model <- survey::svyglm(
80+
#' htn01 ~ sex * insulin + age,
81+
#' design = design, family = quasibinomial()
82+
#' )
10283
#'
103-
#' # --- 3. Calculate all additive interactions ---
104-
#' all_interactions_table <- addintlist(
105-
#' model = interaction_model_logit,
106-
#' factor1_name = "Race1",
107-
#' factor2_name = "ObeseStatus",
108-
#' measures = "all"
109-
#' )
110-
#'
111-
#' # Print the results table
112-
#' print(all_interactions_table, n=50)
113-
#'
114-
#' } else {
115-
#' print("Required packages (survey, NHANES, dplyr, tidyr, msm) not found.")
116-
#' }
84+
#' # RERI, AP and S for every factor-level combination. The Scale column flags
85+
#' # that these are OR-based (see Details for the rare-outcome caveat).
86+
#' addintlist(
87+
#' model = interaction_model,
88+
#' factor1_name = "sex",
89+
#' factor2_name = "insulin",
90+
#' measures = "all"
91+
#' )
11792
#' }
11893
addintlist <- function(model,
11994
factor1_name,
@@ -307,10 +282,29 @@ addintlist <- function(model,
307282
}
308283

309284
# --- Final structure ---
285+
# The Scale column makes explicit which ratio measure RERI/AP/S were built
286+
# from (OR for logistic, HR for Cox, RR for log-binomial/Poisson), since these
287+
# additive-interaction measures are only RR-scale-correct when the outcome is
288+
# rare (see Details).
310289
final_table <- final_table %>%
311-
# Use .data pronoun
290+
dplyr::mutate(Scale = .effect_measure_scale(model)) %>%
312291
dplyr::relocate("Factor1", "Level1", "Factor2", "Level2",
313-
"Measure", "Estimate", "SE", "CI_low", "CI_upp")
292+
"Measure", "Estimate", "SE", "CI_low", "CI_upp", "Scale")
314293

315294
return(final_table)
316295
}
296+
297+
# Internal: infer the ratio scale that RERI/AP/S were computed on, from the
298+
# fitted model's class/family. Not exported.
299+
.effect_measure_scale <- function(model) {
300+
cl <- class(model)[1]
301+
if (cl %in% c("svycoxph", "coxph")) return("HR (Cox)")
302+
fam <- tryCatch(stats::family(model)$family, error = function(e) NA_character_)
303+
link <- tryCatch(stats::family(model)$link, error = function(e) NA_character_)
304+
if (!is.na(fam)) {
305+
if (grepl("binomial", fam) && identical(link, "logit")) return("OR (logistic)")
306+
if (grepl("binomial", fam) && identical(link, "log")) return("RR (log-binomial)")
307+
if (grepl("poisson", fam)) return("RR (Poisson)")
308+
}
309+
"ratio"
310+
}

‎R/jointeffects.R‎

Lines changed: 17 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -47,70 +47,26 @@
4747
#'
4848
#' @examples
4949
#' \donttest{
50-
#' # --- Load required libraries for the example ---
51-
#' if (requireNamespace("survey", quietly = TRUE) &&
52-
#' requireNamespace("NHANES", quietly = TRUE) &&
53-
#' requireNamespace("dplyr", quietly = TRUE) &&
54-
#' requireNamespace("tidyr", quietly = TRUE) &&
55-
#' requireNamespace("msm", quietly = TRUE)) { # msm needed for ratio SE
50+
#' data(nhanes_mortality, package = "svyTable1")
51+
#' nhanes_mortality$htn01 <- as.numeric(nhanes_mortality$htn == "Yes")
5652
#'
57-
#' library(survey)
58-
#' library(NHANES)
59-
#' library(dplyr)
60-
#' library(tidyr)
61-
#' library(msm) # Load msm
53+
#' design <- survey::svydesign(
54+
#' id = ~psu, strata = ~strata, weights = ~survey_weight,
55+
#' nest = TRUE, data = nhanes_mortality
56+
#' )
6257
#'
63-
#' # --- 1. Data Preparation (NHANES Example) ---
64-
#' data(NHANESraw)
58+
#' interaction_model <- survey::svyglm(
59+
#' htn01 ~ sex * insulin + age,
60+
#' design = design, family = quasibinomial()
61+
#' )
6562
#'
66-
#' vars_needed <- c("Age", "Race1", "BPSysAve", "BMI", "ObeseStatus", "Hypertension_130",
67-
#' "SDMVPSU", "SDMVSTRA", "WTMEC2YR")
68-
#'
69-
#' nhanes_adults_processed <- NHANESraw %>%
70-
#' filter(Age >= 20) %>%
71-
#' mutate(
72-
#' ObeseStatus = factor(ifelse(BMI >= 30, "Obese", "Not Obese"),
73-
#' levels = c("Not Obese", "Obese")),
74-
#' Hypertension_130 = factor(ifelse(BPSysAve >= 130, "Yes", "No"),
75-
#' levels = c("No", "Yes")),
76-
#' Race1 = relevel(as.factor(Race1), ref = "White")
77-
#' ) %>%
78-
#' select(all_of(vars_needed)) %>%
79-
#' drop_na()
80-
#'
81-
#' adult_design_binary <- svydesign(
82-
#' id = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~WTMEC2YR,
83-
#' nest = TRUE, data = nhanes_adults_processed
84-
#' )
85-
#'
86-
#' # --- 2. Fit the Interaction Model ---
87-
#' interaction_model_logit <- svyglm(
88-
#' Hypertension_130 ~ Race1 * ObeseStatus + Age,
89-
#' design = adult_design_binary, family = quasibinomial()
90-
#' )
91-
#'
92-
#' # --- 3. Run the function ---
93-
#'
94-
#' # Example 1: Output on Ratio scale (default)
95-
#' joint_effects_ratio <- jointeffects(
96-
#' interaction_model = interaction_model_logit,
97-
#' factor1_name = "Race1",
98-
#' factor2_name = "ObeseStatus",
99-
#' scale = "ratio"
100-
#' )
101-
#' print("--- Joint Effects (Ratio Scale) ---")
102-
#' print(joint_effects_ratio, n = 50)
103-
#'
104-
#' # Example 2: Output on Log scale
105-
#' joint_effects_log <- jointeffects(
106-
#' interaction_model = interaction_model_logit,
107-
#' factor1_name = "Race1",
108-
#' factor2_name = "ObeseStatus",
109-
#' scale = "log"
110-
#' )
111-
#' print("--- Joint Effects (Log Scale) ---")
112-
#' print(joint_effects_log, n = 50)
113-
#' }
63+
#' # Joint effects of sex x insulin on the ratio (OR) scale.
64+
#' jointeffects(
65+
#' interaction_model = interaction_model,
66+
#' factor1_name = "sex",
67+
#' factor2_name = "insulin",
68+
#' scale = "ratio"
69+
#' )
11470
#' }
11571
jointeffects <- function(interaction_model,
11672
factor1_name,

‎R/plotint.R‎

Lines changed: 2 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -74,95 +74,8 @@
7474
#' legend_title = "log(Total Calories)"
7575
#' )
7676
#'
77-
#' # You can now access the plot data if you want it
78-
#' # head(plot_data$plot_data)
79-
#'
80-
#'
81-
#' # =================================================================
82-
#' # Example 2: Continuous x Continuous (NHANESraw)
83-
#' # =================================================================
84-
#' library(NHANES)
85-
#' library(dplyr)
86-
#' library(survey)
87-
#' library(emmeans)
88-
#' library(ggplot2)
89-
#'
90-
#' # --- 3. DATA PREPARATION (NHANESraw) ---
91-
#' data("NHANESraw", package = "NHANES")
92-
#'
93-
#' # Select relevant variables and filter for adults with complete cases
94-
#' nhanes_clean <- NHANESraw %>%
95-
#' filter(Age >= 20) %>%
96-
#' mutate(
97-
#' # Ensure factors are set correctly
98-
#' Gender = factor(Gender),
99-
#' Race1 = factor(Race1)
100-
#' ) %>%
101-
#' # Keep only the variables we need for this model
102-
#' select(
103-
#' BPSysAve, Age, BMI, Gender, Race1,
104-
#' SDMVPSU, SDMVSTRA, WTMEC2YR
105-
#' ) %>%
106-
#' na.omit() # Listwise deletion for simplicity
107-
#'
108-
#' # --- 4. DESIGN & MODEL FITTING ---
109-
#' # Create the survey design object
110-
#' design_nhanes <- svydesign(
111-
#' id = ~SDMVPSU,
112-
#' strata = ~SDMVSTRA,
113-
#' weights = ~WTMEC2YR,
114-
#' nest = TRUE,
115-
#' data = nhanes_clean
116-
#' )
117-
#'
118-
#' # Fit the model: Interaction of Age * BMI on Geometric Mean of BP
119-
#' # We use gaussian(link="log") to model Geometric Means
120-
#' model_bp_interaction <- svyglm(
121-
#' BPSysAve ~ Age * BMI + Gender + Race1,
122-
#' design = design_nhanes,
123-
#' family = gaussian(link = "log")
124-
#' )
125-
#'
126-
#' # --- 5. PLOT THE INTERACTION ---
127-
#'
128-
#' # --- Plot 1: Standard Plot (using quantiles and CIs) ---
129-
#' # This shows the effect of Age (x-axis) at 10th/50th/90th percentiles of BMI
130-
#' print(
131-
#' plotint(
132-
#' model = model_bp_interaction,
133-
#' effect = "Age",
134-
#' moderator = "BMI",
135-
#' data = nhanes_clean,
136-
#' mod_scale = "quantile", # Use 10th/50th/90th percentiles for BMI
137-
#' type = "response", # Plot the Geometric Mean (back-transform from log)
138-
#' show_ci = TRUE,
139-
#' bw = FALSE,
140-
#' pval_position = "bottom.right",
141-
#' xlab = "Age (Years)",
142-
#' ylab = "Predicted Geometric Mean SBP (mmHg)",
143-
#' legend_title = "BMI (Percentiles)"
144-
#' )
145-
#' )
146-
#'
147-
#' # --- Plot 2: Cleaner Plot (using SD and CIs) ---
148-
#' # This shows the effect of Age at Mean +/- 1 SD of BMI
149-
#' # We turn off CIs for a cleaner look if they are too wide
150-
#' print(
151-
#' plotint(
152-
#' model = model_bp_interaction,
153-
#' effect = "Age",
154-
#' moderator = "BMI",
155-
#' data = nhanes_clean,
156-
#' mod_scale = "sd", # Use Mean +/- 1 SD for BMI
157-
#' type = "response",
158-
#' show_ci = FALSE, # Turn off CI ribbons
159-
#' bw = TRUE, # Use B&W linetypes
160-
#' pval_position = "bottom.right",
161-
#' xlab = "Age (Years)",
162-
#' ylab = "Predicted Geometric Mean SBP (mmHg)",
163-
#' legend_title = "BMI (Standard Deviations)"
164-
#' )
165-
#' )
77+
#' # You can now access the underlying plot data:
78+
#' head(plot_data$plot_data)
16679
#' }
16780
plotint <- function(model,
16881
effect,

0 commit comments

Comments
 (0)