Calculates measures of additive interaction (Relative Excess Risk due to Interaction - RERI, Attributable Proportion due to interaction - AP, and Synergy Index - S) with delta method confidence intervals. It can use coefficients from either a joint variable model or an interaction term model (`svycoxph`, `svyglm`, `coxph`, `glm`). It correctly handles the variance-covariance matrix from survey design objects.
Usage
addint(
model,
type = c("joint", "interaction"),
coef_names,
measures = "all",
conf.level = 0.95
)Arguments
- model
A fitted model object (e.g., `svycoxph`, `svyglm`).
- type
Character string: `"joint"` if using a combined categorical variable, `"interaction"` if using main effects and a product term.
- coef_names
A list containing the exact names of the coefficients required for the calculation, based on the model `type`:
If `type = "joint"`: `list(exp1_level = "coef_name_A_only", exp2_level = "coef_name_B_only", both_levels = "coef_name_A_and_B")`
If `type = "interaction"`: `list(exp1_coef = "coef_name_mainA", exp2_coef = "coef_name_mainB", inter_coef = "coef_name_A_x_B")`
Coefficient names must match `names(coef(model))` exactly.
- measures
A character vector specifying which measures to calculate. Options: `"RERI"`, `"AP"`, `"S"`, or `"all"`. Default is `"all"`.
- conf.level
Confidence level for the interval (default 0.95).
Value
A list containing named vectors for each requested measure. Each vector includes Estimate, SE (SE for RERI and AP, SE_log for S), LowerCI, UpperCI. Returns `NULL` or a partial list with `NA` values if calculation fails (e.g., missing coefficients, delta method error).
Details
The function extracts coefficients and the variance-covariance matrix from the fitted model. Based on the `type` specified ('joint' or 'interaction'), it selects the appropriate coefficients corresponding to the exposure levels provided in `coef_names`. It then uses the delta method via the `msm::deltamethod` function to calculate the standard errors for RERI, AP, and log(S), and constructs confidence intervals.
RERI = RR11 - RR10 - RR01 + 1 (or ORs/HRs) AP = RERI / RR11 S = (RR11 - 1) / ((RR10 - 1) + (RR01 - 1))
Confidence intervals for S are calculated on the log scale and then exponentiated.
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)) {
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)
# --- 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.")
}
#> Loading required package: grid
#> Loading required package: Matrix
#> Loading required package: survival
#>
#> Attaching package: 'survey'
#> The following object is masked from 'package:graphics':
#>
#> dotchart
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
#>
#> Attaching package: 'tidyr'
#> The following objects are masked from 'package:Matrix':
#>
#> expand, pack, unpack
#> [1] "--- Results from Joint Model (Black vs Obese Interaction) ---"
#> $RERI
#> RERI_Estimate RERI_SE RERI_CI95_low RERI_CI95_upp
#> 0.02486171 0.31237142 -0.58737502 0.63709844
#>
#> $AP
#> AP_Estimate AP_SE AP_CI95_low AP_CI95_upp
#> 0.009809106 0.122477110 -0.230241619 0.249859830
#>
#> $S
#> S_Estimate S_SE_log S_CI95_low S_CI95_upp
#> 1.0164681 0.2047878 0.6804215 1.5184814
#>
#> [1] "--- Results from Interaction Model (Black vs Obese Interaction) ---"
#> $RERI
#> RERI_Estimate RERI_SE RERI_CI95_low RERI_CI95_upp
#> 0.02486171 0.31237142 -0.58737502 0.63709844
#>
#> $AP
#> AP_Estimate AP_SE AP_CI95_low AP_CI95_upp
#> 0.009809106 0.122477110 -0.230241619 0.249859830
#>
#> $S
#> S_Estimate S_SE_log S_CI95_low S_CI95_upp
#> 1.0164681 0.2047878 0.6804215 1.5184814
#>
# }